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PREFACE 


Contract  F  49620-77-0130  of  the  Air  Force  Office  of 
Scientific  Research,  beginning  September  1,  1977,  was  specifically 
directed  toward  solving  problems  of  elastic  and  viscoelastic  wave 
scattering  and  diffraction  with  applications  to  the  silo  and  its 
related  problems. 

The  Contract  has  been  successfully  executed  at  the  Aldridge 
Laboratory  of  Applied  Geophysics.  This  Final  Report  reflects  the 
accomplishments  of  the  Contract,  including  doctoral  dissertations 
to  be  published  and  unpublished  papers. 

I  am  grateful  to  Mr.  William  Best  of  the  Air  Force  Office 
of  Scientific  Research  for  monitoring  the  contract,  and  Drs.  Ker 
Thomson,  Gerry  Cabaniss,  and  Frank  Crowley  for  their  interest  in 
the  project.  I  wish  to  acknowledge  Ms.  Linda  Ripps,  who  typed  this 
report  most  patiently. 
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GENERAL  INTRODUCTION 

u 

Under  the  Air  Force  Contract,  AF-49620-77-C-0130,  "Elastic 
and  Viscoelastic  Wave  Scattering  and  Diffractions,"  the  main 
objective  of  the  research  was  to  study  the  problem  of  two-  and 
three-dimensional  wave  scattering  and  diffraction  in  elastodynamics 
and  viscoelastodvnamics .  Presently  available  analytical  techniques 
for  solving  wave  propagation  problems  are  useful  only  for  simple 
cases.  In  practice,  the  presence  of  inhomogeneities  and  irregular 
boundary  conditions  defies  analytical  solutions.  One  of  the  best 
numerical  techniques  suitable  for  solving  wave  propagation  in  a 
complex  geological  medium,  such  as  the  problem  of  the  ground  response 
to  seismic  disturbances  in  alluvial  valleys,  is  the  method  of  com¬ 
bining  the  finite  element  method  (FEM)  in  space  and  the  finite 
differences  method  (FDM)  in  time.  The  advantages  of  using  the  finite 
element  method  in  space  are: 

(1)  allowing  for  almost  any  type  of  static,  dynamic,  and  thermal 
loading  to  be  applied. 

(2)  relatively  easy  to  apply  boundary  conditions. 

(3)  its  flexibility  in  modeling  irregular  geology  and  topography. 

(4)  its  distribution  of  errors,  which  are  averaged  over  the  elements 

/( 

throughout  the  domain  in  question. 

The  advantage  of  using  the  explicit  finite  differences  method  in  time 
integration  is  that  solutions  so  obtained  are  generally  conditionally 
stable,  and  fast  convergent  to  yield  solutions  which  are  only  simple 
products  and  vector  additions. 


However,  the  physical  insight  to  the  scattering  and  diffraction 
phenomena  only  can  be  gained  through  analytical  approaches  or 
physical  modeling  studies;  the  interpretation  of  the  numerical  solu¬ 
tions  must  base  on  the  fundamental  understanding  of  the  scattering  and 
diffraction  phenomena.  Therefore,  our  approaches  to  solving  the  com¬ 
plex  scattering  and  diffraction  problems  have  been  oriented  toward; 

(1)  to  solve  practical  problems  by  numerical  methods. 

(2)  to  gain  insight  to  the  numerical  solutions  to  practical  problems 
by  analytical  and  physical  modeling  studies. 

As  for  the  numerical  methods,  we  emphasize  the  research  on  the  finite 
element  formulation  and  computer  programming.  We  have  completed  the 
finite  element  algorithm  for  transient  elastic  and  viscoelastic  wave 
scattering  and  diffraction  problems  in  two  and  thrca  dimensions  and, 
as  a  by-product,  the  finite  element  elastostatic  algorithm  for  static 
and  quasistatic  loading  problems. 

So  far  as  the  effect  of  damping  on  the  elastic  wave  propo- 
gation  is  concerned,  the  computer  codes  for  the  two-  and  three-dimen¬ 
sional  Rayleigh  type  of  damping  have  been  written  and  successfully 
tested.  However,  these  original  computer  codes  have  followed  the 
conventional  "block-by-block"  format  that  require  an  enormously  large 
in-core  storage  and  a  large  amount  of  computing  time  for  solving  geo¬ 
logically  practical  problems. 

Therefore,  we  have  also  been  working  on  the  areas  of: 

(1)  overcome  the  practical  difficulties  of  using  large  in-core  storage. 

(2)  sparse  matrix  without  calculating  zero-matrix  elements  to  save  the 
computing  time. 


iii 


(3)  the  basic  problem  of  convergence  and  stability  to  obtain  a  good 
solution  without  taking  more  computing  time  than  actually  re¬ 
quired. 

(4)  the  undesirable  reflections  from  the  artificial  boundaries. 

The  finite  element  algorithm  for  elastic  waves  has  been 
successfully  applied  to  solve  the  problem  of  the  transient  seismic 
response  of  an  elastic  whole-space,  an  elastic  half-space,  and  a 
layered  elastic  half-space  due  to  a  finite  cylindrical  cavity  source. 

Only  by  a  considerable  reduction  of  in-core  storage  and 
saving  of  computing  time  would  it  be  possible  to  perform  an  extensive 
parameter-s -  ly  of  more  realistic  geological  problems  as  shown  in 
Figure  1. 

In  the  area  of  analytical  studies,  we  have  been  exploring 
one  of  the  most  attractive  analytical  methods,  i.e.  the  high-order 
perturbation  method.  The  existence  of  two  different  types  of  elastic 
waves,  viz.  compressional  and  shear  waves  in  elastodynamic  boundary 
value  problems,  has  frustrated  many  investigators  to  solve  more 
realistic  problems.  In  the  perturbation  method,  we  are  able  to  replace 
the  two  different  wave  numbers  for  P-wave  and  S-wave  by  their  root- 
mean-square  (RMS)  values  so  that  there  will  be  only  one  single  wave 
number  involved.  The  previous  RMS  perturbation  solutions  developed  by 
Pao  and  Thau  (1967)  are  only  up  to  the  order  of  two.  We  have  now 
extended  the  method  to  the  higher-order  solutions,  through  the  general 
recurrence  formulas.  Moreover,  we  have  applied  the  newly  developed 
perturbation  method  to  solve  one  of  the  long-standing  elastodynamic 
wedge  problems  and  obtained  very  encouraging  results.  In  order  to 
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MORE  REALISTIC  GEOLOGICAL  STRUCTURES 


gain  further  insight  to  wave  scattering  and  diffraction  problems, 
three-dimensional  physical  model  experiments  have  been  performed 
at  the  Aldridge  Laboratory  of  Applied  Geophysics  to  study  the 
acoustic-elastic  wave  scattering  and  diffraction  by  a  vertical 
elastic  cylinder  immersed  in  a  fluid  due  to  a  transient  acoustic 
point  source. 

Physical  insight  into  these  problems,  both  analytically 
and  experimentally,  has  been  invaluable  in  interpretation  of  the 
finite  element  results. 


ABSTRACT 


The  three-dimensional  finite  element  for  both  elastic  and 
viscoelastic  wave  scattering  and  diffraction,  based  on  the  principle 
of  virtual  work, has  now  been  completed  for  solving  transient  elasto- 
and  viscoelasto-dynamic  problems,  as  well  as  elasto-  and  viscoelasto- 
static  problems.  Computer  programs  for  the  finite  element  formula¬ 
tion  are  developed  and  tested. 

The  conventional  "block-by-block"  format  used  in  the  computer 
programs  has  been  improved  in  order  to  save  the  in-core  storage  and 
computer  time.  Stability  condition  and  convergence  criteria  of  the 
finite  element  in  space  and  the  finite  differences  in  time  have  been 
examined.  The  computer  programs  developed  under  the  Contract  are 
now  ready  and  have  been  submitted  to  AFGL  for  their  use. 
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FINITE  ELEMENT  FORMULATION 


Based  on  the  principle  of  virtual  work,  we  have  now 
completed  the  three-dimensional  finite  element  formulation  for 
viscoelastic  wave  scattering  and  diffraction  to  complement  that 
for  elastic  wave  scattering  and  diffraction  due  to  a  transient 
seismic  disturbance,  as  reported  in  the  Final  Report  of  an  earlier 
contract  AFOSR-76-2968 ,  "Numerical  and  Analytical  Solutions  to 
Elastodynamic  Problems".  The  general  procedure  for  using  the  finite 
element  method  in  solving  elastodynamic  problems  in  continuum  mechanics 
includes  idealization  of  the  problem,  space  discretization  of  finite 
elements,  and  formulation  of  the  set  of  simultaneous  equations  to  be 
solved.  For  viscoelastic  waves,  the  equations  of  motion  for  elasto- 
dynamics  must  include  a  term  to  account  for  energy  dissipation. 

Usually  this  mechanism  takes  the  form  of  viscous  damping,  which  is 
linearly  proportional  to  velocity.  Hence,  the  equations  of  motion 
for  viscoelastodynamics  are: 

(1)  [m]{u}  +  [c]  {u}  +  [k]  lu}  =  { f (t) } 

where  {u},  {u},  {u}  are  the  displacement,  velocity  and  acceleration 
vectors  of  the  finite  element  assemblage,  and 
[m]  =  the  assembled  mass  matrix. 

[c]  =  the  viscous  damping  matrix. 

[k]  *  the  assembled  structural  stiffness  matrix. 

{f(t)}  ■  the  external  time-dependent  vector  load  vector. 
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It  is  obvious  that  equation  (1)  reduces  to  the  equations 
of  motion  for  elastodynamics  as  [c]  =  0.  Equation  (1)  states  that, 
in  principle,  static  equilibrium  at  time  t,  which  includes  the 
effects  of  inertia  force,  elastic  force  and  damping  force,  is 
reached.  While  in  static  analysis,  one  simply  neglects  the  inertia 
and  damping  effects  in  equation  (1). 

The  free  vibration  equilibrium  equation  for  the  undamped 
system  is  well-known  as  follows: 


(2) 


[m]  {u}  +  [k]  {u}  =  0 


Assume  that 


(3)  {U}  =  {$}  T (t) 


where 


(4)  T(t)  =  A  sin  a  t  +  B  cos  oj  t 

{ }  =  time  dependent  shape  vector,  and 
A  and  B  =>  arbitrary  constants.  Differentiating  equation  (3) 
twice  with  respect  to  time  t,  we  have 


(5) 


{u}  =  {<J>}  T (t)  =  -  a)2  {$}  T(t) 


.  ..  - 
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Equation  (2)  yields 

(6)  ([k]  -  w2  [m] )  { <J»>  =  0 

In  order  to  have  a  non-trivial  solution  to  this  equation,  the 
following  condition  must  be  satisfied 

(7)  [  [k]  -  u2[m]  j  =  p(uj2)  =  0 

This  leads  to  the  classical  eigenvalue  problem.  The  solution  to 

this  equation  is  in=tD^,  n=l,2,3 . .  N  in  which  N  is  the  number 

of  the  modes,  are  the  natural  angular  frequencies  of  the 
system.  The  shape  vector  function  {>}>}  can  be  obtained  by  solving 
equation  (6)  for  each  value  of  ui,  giving 

(8)  [k]  {<j>n}  =  u2  [m]  { <pn> 

for  the  n-th  mode  of  vibration. 

Superposition  of  the  response  for  all  modes  gives 

N 

(9)  {U}  =  E  U}{T(t)}  =  [*]  {T (t)  } 

n=l  n  n 

where  the  n-th  column  in[$]  contains  the  shape  vector  {$  } 
Substituting  equation  (9)  into  (1) ,  we  have 


(10)  [m]  [*]  {T}  +  [c]  [*]  {T}  +  [kj  It]  {T}  =  {f (t) } 


Premultiplying  equation  (10)  by  f^]  ,  we  have 

(11)  [$]T  [m]  [*]  {T}  +  [<(>]T  [c]  [*]  {T}  +  [J>]T  [k]  [$] 

=  [$]T{f(t)}=  (F(t) } 

The  Orthogonality  conditions  are: 


{V  [m]  {$n}  =  0 


(12) 


if  m  ?  n 


Umr  [k]  {$n}  =  0 


Caughey  (1960)  pointed  out  that: 

T 

a)  Equation  (11)  will  be  decoupled  if,  and  only  if,  [<*>]  [c]  [4>]  is 

diagonal, 

T 

b)  if  [$]  [c]  [<t]  is  diagonal,  the  damped  system  will  have  normal 

modes,  which  are  identical  with  those  for  the  undamped  system,  and 

c)  if  a  damping  matrix  [C]=[$]  tc] [<J>]  has  the  form 


(13) 


N  -1  ,  b 

[C]  =  [M]  Z  a.  ([M]  1  [K])° 

fcRO 


the  damped  system  will  have  classical  normal  modes. 
Now,  if  the  damping  matrix  also  satisfies 


(14)  {$m)  [C]  {<f>n}  =0  for  m/n 


and  equation  (11)  becomes 


(15) 


M. 


L  T  +  C  T  +  K  T  =  a  (t) 
n  n  n  n  n  n  ti 


t 
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where 


Mn  =  Un>T  [ml  { <Pn> 
(16)  Cn  =  {^n}T  lc]  Un} 

Kn  -  (,/  !*](*„) 


Since  Mn,  Cn>  Kn  are  scalars,  equation  (1)  is  decoupled  to  N  linear 
equations.  Each  of  these  equations  represents  the  equation  of  motion 
of  a  single-degree-of-freedom  system.  The  damping  constant  can  be 
so  chosen  that 


(17) 


C  =  2  X  w  rn 
n  n  n  n 


(18) 


“n-  'W 


where  X^=  the  damping  ratio  for  the  n-th  mode. 

«n=  the  natural  angular  frequency  for  the  n-th  mode.  Then,  the 
damping  ratio  in  each  mode,  X^  can  be  related  to  the  natural  angular 
frequency  of  that  mode  by 


N-l 


2b-l 


(19)  A  =  -  L  a.  w_ 

n  2  b=0  b  n 


In  equation  (13),  if  N=2, 


(20)  [C]  =  ao  (M)  +  ax  [K] 
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This  form  was  first  proposed  by  Rayleigh.  Equation  (19)  becomes 


(21) 


A  = 
n 


+  al 


Therefore,  Rayleigh  damping  can  be  considered  as  a  special  case  of 
Caughey  damping.  The  constants  aQ  and  a^  can  be  determined  by  any 
two  arbitrary  conditions.  If  these  two  conditions  are  so  chosen 
that  the  damping  ratios  are  specified  for  the  first  two  modes,  then 


2((i)-A1  -  u) .  A « ) 

(22)  a  =  - — - — - 

O  0)  -  <D 

(  -  m2  ) 


(23) 


2  (  A2/co1  -  XVw2) 

OJ  X  -1 

(  2/o)1  -  1/w2  ) 


The  procedure  described  above  involves  the  calculation  of  the 
frequencies  for  the  first  two  modes. 

Another  procedure  for  determining  the  values  of  and  a^ 

is  to  specify  the  damping  ratio  at  the  frequency  for  which  the  damping 
is  a  minimum,  usually  the  fundamental  frequency  of  the  system.  Then 
we  have 


(24) 


ao  *  Ai 


(25) 


a1  *  Al/ul 
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4 


> 
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Solving  the  eigenvalue  of  the  equation 

(7)  |  [k]  -  u)2  [m]  |  =  p(oj2  )  =0 

is  equvalent  to  calculating  the  roots  of  the  polynomial  p(ut^),  which 
has  the  order  equal  to  the  order  of  [k]  and  [m]  that  are  usually  a 
very  large  number.  Therefore,  an  iteration  solution  method  has  to  be 
used.  The  time  integration  scheme  with  damping  effect  then  is 


(26)  {u (t+At) }  =  {u(t)}  +  {u(t)}At 

{u(t+At)}  =  {u (t) }  -  [M]  1  [K]  {u (t+At)} At 

-  [M]_1[C]{u(t) }At  + [M] _1  {F (t) }  At 


(27) 


TEST 


Based  on  the  foregoing-  formulation,  we  have  developed  the 
computer  programs.  As  an  initial  test  of  our  three-dimensional  program, 
we  have  studied  the  problem  of  wave  propagation  in  a  viscoelastic 
rectangular  bar  with  one  end  of  the  bar  fixed.  The  free  end  of  the 
bar  is  subjected  to  a  uniform  external  force  load.  The  simulated 
forcing  function  is  that  of  the  sinusoidal  type,  using  a  four  time- 
step  of  increasing  (or  decreasing)  time.  The  maximum  of  the  magnitude 
for  the  forcing  function  for  the  test  case  is  lO'*'0  dynes.  In  this 
test  problem,  we  used  only  20  hexahedron  elements.  Use  of  symmetry 
is  made  to  reduce  the  size  of  the  problem  and  the  boundary  displace¬ 
ments  are  prescribed  in  a  manner  shown  in  Figure  2.  The  elastic 
constants  used  are:  Poisson's  ratio  =  0.25,  Young's  Modulus  =  4. 5x10 10 
dynes/cm  .  The  dynamic  damping  behavior  of  unconsolidated  earth 
materials  primarily  depend  on  the  parameter  of  damping  ratio  Xn, 
which  must  be  determined  by  in-situ  experimental  measurements.  On 
the  basis  of  spectral  studies  of  the  damping  behavior  of  unconsolidated 
earth  materials,  the  U.S.  Atomic  Energy  Commission  Regulatory  Guide, 
Directorate  of  Regulatory  Standards  indicates  that  the  reasonable 
values  of  parameter  \n  are  approximately  in  the  range  of  0.005  to  0.1. 
Figures  3  and  4  show  the  displacement-component  responses  on  the 
surface  of  the  bar,  the  observation  point  is  located  at  the  nodal  point 
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of  71  as  shown  in  Figure  2,  for  the  three  cases: 


(i)  Response  without  damping. 

(ii)  Response  with  a  damping  ratio  A^  =  0.02. 

(iii)  Response  with  a  damping  ratio  A^  =  0.05. 

The  responses  of  and  are  exactly  the  same  due  to  symmetry. 

In  Figures  3  and  4,  the  magnitudes  of  Ux  (or  Uy)  and  Uz  are  drastically 
reduced  for  both  the  cases  of  A  =  0.05  and  =  0.02.  In  this  particu¬ 
lar  example,  there  is  a  slight  phase  change  in  Uv  (or  U„)  component, 

x  y 

while  the  phase  change  in  the  Uz  component  is  apparent. 


ions  of  an  Elastic  Bar  and 


MODIFYING  ELASTODYSTATIC  AND  ELASTODYNAMIC 
FINITE  ELEMENT  COMPUTER  PROGRAMS  FOR 
REDUCING  THE  IN-CORE  STORAGE 

As  we  mentioned  before,  our  original  finite  element  computer 
programs  have  been  developed  in  the  conventional  "block-by-block” 
format  so  that  a  large  in-core  storage  for  the  global  stiffness 
matrix  was  required  and  a  large  amount  of  computing  time  was  also  re¬ 
quired.  If  we  use  the  "block-by-block"  format,  it  is  a  formidable 
task  to  handle  even  a  small  practical  problem  of  a  more  realistic 
geological  structure  by  the  present  art  of  the  largest  digital 
computers  now  available. 

In  the  past  year,  we  have  substantially  modified  our  three- 
dimensional  elastoaynamic  and  elastostatic  computer  programs  for 
solving  more  realistic  geological  problems.  We  have  overcome  many 
serious  difficulties  of  the  finite  element  method  for  both  elasto- 
dynamic  and  elastostatic  problems. 

(A)  Elastodynamics 

The  final  form  of  the  finite  elemerb  forum.1  I’-ion  of  the 
equations  of  motion  in  elastodynamics  is 

(28)  [Ml {u}  +  [K]  {u}  =  (F(t) } 

A  Gaussian  elimination  algorithm  is  generally  used  to  solve  a 
discrete  system  of  equations.  However,  matrices  [K]  and  [M]  must 
be  in-core  stored  in  the  calculation  that  creates  stringent  limitation 
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on  the  total  number  of  equations  to  be  solved. 

Because  the  matrix  [K]  is  usually  sparse,  a  special  ordering 
of  the  nodal  points  can  produce  a  matrix  [K] ,  which  is  tightly 
banded  about  the  main  diagonal.  In  our  original  version  of  the 
Aldridge  Finite  Element  Algorithm  (AFEA) ,  we  used  the  conventional 
"block-by-block"  format  with  the  number  of  equations  per  block 
known,  the  stiffness  matrix  is  assembled  in  one  block  at  a  time  and 
is  stored  on  the  scratch  disk.  This  "block-by-block"  computer 
program  requires  (2B,  B)  words  in-core  for  [K]  ,  where  B  is  the 
half-band  width.  For  example,  for  a  three-dimensional  50x50x50 
element  structure,  we  would  have  to  solve  some  397,953  simultaneous 
equations  for  the  assembly  stiffness  matrix  with  a  half-band  width 
of  B=7959 ,  and  an  in-core  storage  of  2x7959x7959x4=506,766  kilobytes. 
The  banded  and  blocked  character  of  [K]  is  as  shown  in  Figure  5. 

The  assembly  mass  matrix  of  the  structure  is  formulated  by  lumping 
the  weighted  element  masses  at  the  nodal  points  of  the  elements. 

The  mass  of  an  element  is  weighted-distributed  to  the  nodal  points. 

In  the  present  modified  version  of  AFEA,  taking  the  advantage  of 
using  the  lumped  mass  matrix,  and  the  finite  differences  time  inte¬ 
gration,  we  thus  avoid  the  conventional  Gaussian  elimination 
algorithm.  Instead,  we  solve  the  following  equations, 

(29)  {u(t+At))N  =  {u(t)lN  +  ({u(t)}At)N 

(30)  {u(t+At)}N  =  {u(t))N  +  (  [M]_1[K){u(t4At)}  )N 

+  (  [M]"1  {F  (t) }  At)N 
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where  N  =  the  total  number  of  equations. 

We  operate  only  on  the  product  {  [M]  ^  [K]  {u(t  +  it)}  )  . 

N 

Thus,  we  need  only  "N"  words  in-core,  no  matter  how  large  the  half¬ 
band  width  B  may  be.  Figure  6  shows  a  comparison  of  the  in-core 
storage  for  the  global  stiffness  matrix  required  between  the  original 
and  present  modified  versions  cf  AFEA . 

Additionally,  the  present  modified  AFEA  codes  save  accounting 
time,  since  only  the  non-zero  matrix  elements  in  [K]  are  calculated. 

The  blocked  matrix  [K]  is  sparse.  The  larger  is  the  half-band  width 
B,  the  more  sparse  the  blocked  [K]  matrix  will  be.  In  the  three-dimen¬ 
sional  case,  the  maximum  non-zero  elements  in  each  row  are  81  for 
hexahedron  elements.  In  two-dimensional  case,  the  maximum  non-zero 
elements  in  each  row  are  18  for  quadrilateral  elements.  In  other 
words,  there  are  (7959-81)  zero  matrix  elements  in  each  row  for  the 
50x50x50  element-structure.  Therefore,  in  using  the  present  modified 
AFEA,  the  computing  time  can  be  reduced  for  solving  larger  element- 
structures  without  the  increasing  I/O  time.  Table  I  gives  the  com¬ 
parison  of  the  accounting  time  for  the  original  and  present  modified 
versions  of  AFEA  in  solving  both  two-dimensional  and  three-dimensional 
elastodynamic  problems  on  an  IBM  360/91. 

As  our  first  example,  in  using  the  modified  three-dimensional 
elastodynamic  program  of  AFEA,  we  studied  the  classic  Lamb's  problem, 
i.e.,  the  problem  of  a  homogeneous,  isotropic,  elastic  half-space, 
subjected  to  an  external  vertical  force  on  the  free  surface,  which 
is  assumed  to  be  a  bell -shaped  function.  In  this  test  problem,  we 
used  12  hexahedron  elements  in  each  dimension  (Figure  7) .  Use  of 
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Versions  of 
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nd  Fix  Boundary 


9.  Comparison  of  the  Finite  clement  and  the  Analytical 
Solutions.  (Finite  Clement  Mesh  :  12x12x12  ) 


symmetry  is  made  to  reduce  the  size  of  the  problem.  The  symmetric 
boundary  conditions  in  Figure  7  are: 


U  =  0  on  OBDC 
x 


U  =  0  on  OBFE 

y 


In  addition  we  used  Smith's  technique  (1974)  to  eliminate  the  single 
boundary  reflections  by  superimposing  the  solutions  of  free-artifi¬ 
cial-boundaries  and  rigid-artificial-boundaries  problems.  Although 
the  technique  fails,  when  multi-reflections  occur  at  the  same 
boundary,  and  has  the  disadvantage  of  requiring  2°  solutions  for  the 
n  number  of  non-reflecting  boundaries,  we  prefer  this  method  at 
present  because  it  offers  the  advantages  of  perfect  elimination,  i.e. 
independent  of  both  incident  angle  and  frequency.  The  Smith's 
technique  is  particularly  advantageous  for  three-dimensional  finite 
element  problems  since  the  number  of  elements  required  is  reduced 
by  8.  Figure  8  shows  the  superposition  of  the  solutions  of  the 
free  and  boundary  problems.  Figure  9  shows  the  comparison  of  the 
finite  element  and  the  analytical  solutions.  The  test  results  are 
highly  accurate  with  the  displacement  component  being  within  about 
2  %  of  the  analytical  value. 


(B)  Elastostatic  Case 

For  the  elastostatic  case,  a  conventional  Gaussian  elimination 
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is  needed.  In  the  new  version  of  AFEA,  we  developed  a  small  band 
solver  that  can  be  used  with  tape-storage  to  solve  a  large  system 
of  linear  equations.  Instead  of  requiring  the  in-core  storage  of 
(2B,  B)  words  for  [K]  in  the  orginal  version,  the  required  in-core 
storage  now  for  [K]  reduces  to  (2  x  NF,  B)  words  where  NF  is  the 
numbers  of  degree  freedom,  that  is,  an  in-core  storage  of  (6,  B) 
words  for  the  three-dimensional  case  and  of  (4,  B)  words  for  the 
two-dimensional  case.  Figure  10  gives  a  comparison  of  the  in-core 
storage  for  the  global  stiffness  matrix  required  in  the  original 
and  new  modified  versions.  It  is  expected  that  the  large  number 
of  calls  to  tape  requires  a  large  amount  of  I/O  time. 


B  (Half-band 


Figure  10.  Comparison  of 
Versions  for  In-core 
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STABILITY  CONDITION  AND  CONVERGENCE  CRITERIA 

OF  THE  FINITE  ELEMENT  IN  SPACE  AND 
THE  FINITE  DIFFERENCES  IN  TIME 


Ever  since  the  finite  element  method  emerged  in  various 
branches  of  engineering  and  physics,  many  investigators  have 
attempted  to  establish  the  convergence  criteria  and  the  stability 
condition  techniques  for  the  FEM  in  space  and  the  FDM  in  time. 

Based  on  the  principle  of  virtual  work,  the  FEM  involves 
the  whole  region  of  interest,  which  is  discretized  into  a  finite 
number  of  subregions.  The  accuracy  of  the  analysis  depends  mainly 
on  the  number  of  elements  used,  and  on  the  nature  of  the  assumed 
displacement  functions  with  the  elements.  The  accuracy  of  analysis 
can  be  increased  by  using  more  elements  in  the  representation  of 
the  structure  provided  that  the  elements  satisfy  certain  conver¬ 
gence  criteria,  that  is,  the  elements  must  be  complete  and  com¬ 
patible. 

Whether  a  specific  element  is  complete  and  compatible 
depends  on  the  formulation  that  must  be  studied  individually.  The 
aim  in  the  numerical  integration  of  the  finite  element  system  equi¬ 
librium  is  to  evaluate  a  good  approximation  to  the  actual  dynamic 
response  of  the  structure  under  consideration.  The  choosing  time 
step  must  correspond  to  the  smallest  period  in  the  system  which 
could  mean  very  small  time  steps. 

In  the  AFEA ,  i.e.,  finite  element  in  space  and  finite  differ¬ 
ences  in  time,  we  found  that  the  elements  of  4CST  (four  constant 
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strain  triangles)  for  two-dimensional  cases  and  5CST  (five  constant 
strain  tetrahedra)  for  three-dimensional  cases  in  space  integration 
are  complete  and  compatible,  and  central  differences  of  time  steps 
in  time  integration  are  adequate.  It  follows  that  the  element  size 
Ax,  Ay,  or  Az  in  comparison  with  the  wavelength  of  wave  propagation 
plays  a  crucial  role  in  convergence.  As  in  time  integration,  the 
number  of  operations  required  is  directly  proportional  to  the  number 
of  time  steps  required  for  the  solution.  The  selection  of  an  appro¬ 
priate  time  step  At  in  time  integration  is  of  great  importance  in 
obtaining  a  stable  solution.  On  one  hand,  the  finite  element  size 
and  the  time  step  must  be  small  enough  to  obtain  a  stable  solution, 
and,  on  the  other  hand,  the  element  size  must  not  be  smaller  than 
necessary,  because  this  would  mean  that  the  computation  is  more 
costly  and  requires  more  in-core  storage  than  actually  needed. 

Choosing  a  time  step  smaller  than  necessary  not  only  increases  the 
cost  of  calculation;  it  also  causes  an  inaccurate  result.  A  common 
mistake  is  choosing  a  coarse  grid  of  elements  and  a  small  time 
step.  We  have  tested  the  convergence  criteria  and  stability  condition 
by  using  one  of  our  computer  program  codes  for  the  two-dimensional  plane 
strain  elastodynamic  calculation,  again  the  best  case  is  that  of 
Lamb's  problem.  The  following  forcing  function  is  used, 

F  (t)  =  ^r-  [H(t-At  )  -  H(t-  3  At)] 

(31)  +  — rrr-[  At2  -  2(t-At  )2]  [H(t-3  At)  -  H(t  -  4  At)] 

uT 

1  2 

+  ~t'2~  [(t  -  4  At)  ]  [H  (t  -  4  At+  At)] 
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in  which  the  factor  At  controls  the  rise  of  the  step.  We  found 

that  with  various  element  sizes  (Ax,  Ay) ,  various  time  steps  (At) , 

and  various  forcing  function  sharpness  factors  (At),  the  present 

finite  element  model  requires  the  ratio  £  =  V^At/As  be  greater 

than  1.2  and  At  be  greater  than  10  At,  where  is  the  com- 

pressional  velocity  and  As  is  the  smallest  width  of  the  element  in 

the  finite  element  assemblage.  As  an  estimate  of  As  is  required, 

it  appears  that  if  the  average  period  is  T  ,  As  would  have  to  be 

about  T  V  /10.  Figure  11  gives  the  finite  element  mesh  of  the  two 
n  p 

dimensional  Lamb's  problem.  In  Figure  11,  use  of  symmetry  is  made 
to  reduce  the  size  of  the  problem.  Although  the  grid  sizes  used 
are  regular,  both  element  size  and  shape  may  vary.  Synthetic 
seismograms  at  various  nodal  points  along  the  free  surface  have 
been  generated  with  5  *  1.5  and  AT  =  12At  (Figure  12a  and  12b). 
Figure  13  shows  the  diagram  of  the  particle  motion  at  an  observation 
point  on  the  free  surface.  The  characteristic  elliptic  retrograde 
particle  motion  of  Rayleigh  wave  is  clearly  identified. 
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II.  THE  PERTURBATION  METHOD  APPLIED  TO  SOLVE 


ELASTODYNAMIC  PROBLEMS 


ABSTRACT 


The  root-mean-square  perturbation  method  for  solving  wave 
scattering  and  diffraction  introduced  by  Thau  and  Pao  has  been 
further  developed  to  include  the  higher  order  terms.  The  method  is 
particularly  suited  for  solving,  for  example,  the  classic  problem 
of  wave  propagation  in  a  wedged  medium.  As  an  illustration,  the 
formal  solutions  to  the  problem  of  SH-wave  propagation  in  a 
wedged  half-spce  are  obtained.  Future  numerical  evaluation  of  the 
solutions  is  in  order. 
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INTRODUCTION 


In  solving  elastodynamic  boundary  value  problems,  the 
dual  existence  of  two  different  types  of  elastic  waves,  compres- 
sional  wave  and  shear  wave,  even  in  a  homogeneous,  isotropic, 
elastic  solid  body  has  frustrated  many  theoreticians  to  solve  some 
practical,  as  yet  complex  problems  in  seismology.  For  example, 
the  use  of  Kontorowich-Lebedev  transform  will  yield  a  solution  to 
the  single  scalar  field  problems  of  acoustic,  SH-  or  electromag¬ 
netic  wave  propagation.  However,  the  use  of  Kontorowich-Lebedev 
transform  is  no  value  in  obtaining  formal  solutions  in  the  elastic 
solid  wedge  problem  of  P  and  S  wave  propagation. 

Thau  and  Pao  (1967)  introduced  a  perturbation  method,  in 
which  the  two  different  wave  numbers  in  the  two  steady  state  wave 
equations  are  replaced  by  their  root-mean- square  (RMS)  average 
value.  All  the  perturbed  wave  equations  and  boundary  conditions 
involve  only  the  RMS  value,  that  is,  there  is  only  one  single 
wave  number  involved.  The  perturbation  parameter,  which  is  a  small 
quantity,  is  the  difference  of  the  actual  wave  number  and  the  RMS- 
value.  They  generated  the  perturbation  solutions  of  the  first 
three  orders,  viz.  the  zeroth,  the  first,  and  the  second  order.  They 
also  illustrated  the  application  of  this  perturbation  method  to  solve 
elastic  wave  diffraction  problems.  Based  on  their  two-term  pertur¬ 
bation  solutions  of  the  zeroth,  and  the  first  order,  the  results  are 
in  good  agreement  with  the  known  exact  near  field  solutions  at  low 
frequencies . 
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In  this  paper,  we  present  the  recursive  formulas,  which 
can  generate  any  higher  order  solutions  for  this  RMS-value  pertur¬ 
bation  method. 


Thau  and  Pao  showed  that,  for  the  steady  state  motion,  the 
two  Helmholtz  equations  for  elastic  waves  are 


(  V2  +  k2 


(  V2  ♦  k2 


where  is  the  Laplacian  operater. 


k1  =  ui  /  v± 


k_  =  a)  /  v. 


w  is  the  frequency,  v^  and  v 2  are  velocities  and  v  >v^ .  For  the 
SH-wave  problems  in  two  media,  v^  and  v 2  are  the  shear  velocities, 
and  w^  and  w^  are  the  z-direction  displacement  components  for  the 
two  different  media,  respectively.  For  waves  in  an  elastic  solid, 
v^  is  the  shear  wave  velocity,  v^  is  the  compressional  wave  velocity, 
and  w^  and  are  the  displacement  potentials.  For  simplicity  and 
without  loss  of  generality,  here  only  one  component  of  the  vector 
displacement  potentials  is  considered.  For  example,  for  the  plane- 
strain  problems,  the  displacement  vector  is 


u  =  ^w2+  £Xkwl 
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(Uniting  the  factor  exp(-io)t),  we  can  easily  find  the  equivalencies 
for  the  acoustic  and  electromagnetic  waves. 

Let 

(6)  k2  =  (kj2  +  k22)/2 

(7)  k  2  =  k2(l  +  2e) 

(8)  k22  =  k2(l  -  2c) 


where 


(9) 


e  = 


2  2 
kl  '  k2 


2  (k2  +  k2) 


is  the  perturbation  parameter 


Assume  that 


(10)  w,  = 


w{o)+  e  w<1J  +  e2w<2)  + 


n  (  n  )  , 

+  e  w  +  ■ 


(11)  w2  -  w*0l»  e  “211  *  c2"22>  * 


Substituting  (7)  ,  (8)  ,  (10)  ,  and  (11)  into  (1)  and  (2)  and 
equating  coefficients  of  like  powers  of  the  n-th  order  wave  equations, 
we  heves 


(12) 


(  V2  +  k2  >w{n) 


,  .2  (n-1) 

-2k 


S 


113) 


(  V2  *  K2  ,.‘nl  -  *  2  k2  v'-1' 


with  n  =  1,2,3, . 

For  n  =  0,  the  wave  equations  become: 


(14) 


(  V2  +  k2  )wJo) 


=  0 


(15) 


(  V2  +  k2  )w<°> 


=  0 


For  higher  orders  (n  /  0) ,  the  total  solutions  of  (12)  and  (13) 
include  two  parts: 


(16) 


(n) 

'lc 


+ 


w 


(n) 

IP 


(17) 


( n 
2 


c 


+ 


w 


( n ) 
2p 


where  w^"^ and  w^^  are  the  complementary  solutions,  which  satisfy  the 
homogeneous  equations 


(18)  (  V2  +  k2  )w<n)  =  0 

lc 

(19)  (  V2  +  k2  =  0 
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w,(n)  and  are  the  particular  solutions,  which  can  be  obtained 

lp  2p 

from  the  lower  order  solutions  and  satisfy  the  following  inhomo¬ 
geneous  equations: 


.  2  (n) 

(V  +  k  )  w. 

IP 


2  k2  w}"'1* 


,  „2  ^  .  2  ,  (n) 

(7  +  k  )  w 

2p 


-  .2  (n-1) 

2  k  w2 


The  particular  solutions  for  the  first,  the  second,  and  the  third 
orders  are: 


r  •  Vw 

a.  a. 


Vw  (o) 

n.  1 


-  r  •  Vw. 


(2)  _  (1)  (1)  1.22  (o) 

(24)  w,  =  r  •  Vw,  -  w,  +  -■  k  r  w, 

IP  0/  A/  1  i-P  2  1 


(25)  w' 


_  (1)  A(ll  1  .2  2  (o) 

’  lW  2  +  w2p  +  ~r  k  r  w2 


{r-  V  [2w.|2^  +wj2)  +w{ 1}  ]  ”4w  ^  +  %  kVw'1'} 

3  %  'v,  1  lc  lc  lp  1 


-{r  •V[2w‘2,+w<2,-w<1)]W2)  +  %  kVw^} 

3  ^  ^  2  2c  2c  2p  2 


The  n-th  order  particular  solutions  are: 


(28)  w 


where 


t  „  ,  r.  ,,  (n-1)  (n-1)  .  (n-2) , 

-  {  r  .7  )  [  (n-1) w  +  aw  +  J 

a^  %  %  1  2  lc  2  lc 

r ,  _  (n-1)  ..22  (n-2) 

-[(2n-3)+b2J  w^  +  *i  k  r  w1 


+  (-1)  [IL  +  I2  +• 


+  I  .  +  I  ]} 
n-1  n 


(29)  w. 


1  r  ,  „  .  r  .  (n-1).  (n-1)  (n-2), 

—  {-(  r  -V  )  [(n-l)w2  +  a2w2c  -  ] 

+  [  (2n-3) +b2]  w2^  1)+  >j  k2r2w2n  2' 


+  [I1  -  Z2  + 


+  (-1) n_1I  .  +  (-l)ni  ]} 
n-1  n 


a^  =  £  (m)  ,  a2  =  £  (m)  ,  a 3  =  £  (m)  ,  a4=  £  (m) 


b  =  £  ( 2m-l ) ,  b  =  £  (2m-l) ,  b  =  £  (2m-l) 

m=l  m=l  m= 1 


I,  =  (r.  V)  (r.V-2)  [a  w1(n_2) +b  w(n-3)  ]  /2  ! 

X  >\,  0/  'X/  %  o  1C  J  1C 


I,  =  (  r  •  V  )  (  r  •  V  -2)  (r-7-4)  [a  .W?(n"3) +b.w.(n'4)  ]/3  ! 
2  *\j  %  ty*  'X,  %  4  lc  4  lc 


(34)  I 


[r . V-2 (n-4) ] 


n- 1 


(r -V)  (r .V-2)  (r-V-4) 

x,  %  a,  ^  a,  *\, 


'x. 


X  [3  W,(3)  +  4w.(2)  ]  /  (  n-  3  )  ! 

1  c  X  c 


t 


i 


(35)  I  =  (r-V) (r.V-2) (r-V-4) 

^  'b  %  %  %  *X*  'X/ 


[r . V-2 (n-4) ]  x 

%  ru 


X 


(2)  A 
w .  + 
lc 


w(I)] 
wlc  J 


/  (  n-  2  )  ! 


1. 


c 


In  obtaining  the  recurrence  formulas  (28)  and  (29) ,  we  used  the 
following  indentity: 


(36)  ( 7  2  +  k  2 )  [  (r.V)  (r.V-2)  (r-V-4) 

0/  \  a.  't  'Xj  a. 


(r.V-2m)  wfn)] 
a,  a,  1C 


-  2 (m+1) k  (r-V) (r-V-2) 


[r.V-2 (m-1) ]  wfn) 

•v  %  1C 


where  m  =  0,1, 2, 3, . .  i  =  1,2. 

It  is  interesting  to  point  out  that  the  first  order 
particular  solutions  presented  here  are  the  same  as  those  which 
were  obtained  by  Thau  and  Pao.  While,  the  second  order  solutions 
are  different  from  the  solutions  obtained  by  Thau  and  Pao.  In  fact, 
we  can  easily  find  three  different  forms  of  the  perturbation  solu¬ 
tions  for  the  second  order  particular  solutions,  viz.. 


(37) 

v;2)= 

ip 

-  r.vt«<0)- 

'V  'X* 

w(1>] 
lc  J 

-  h 

(38) 

<2) 

w_  = 

-  r*7[w^0^  + 

w!15] 
2c  J 

-  h 

2p 

a,  a.  2 

.2  2  (o) 

k  r  w^ 


.2  2  (o) 

k  r  w2 


9 


t 


c 


(39) 


(2)  n  ,  (1)  (o)  ,  ^  ,  .2  2  (o) 

=  r  •  V  [  w ,  -  w ,  J  +  h  k  r  w. 

Ip  ^  'v  1  1  1 


(40)  w{2)  =  r-V  [wi15  +  w*o)]  +  >j  k2r2  v/o) 

2p  ^  2  2  2 


(41) 


(2) 

'IP 


r  •  7  [Wl(1) 
%  lc 


(42) 


i[2) -  -  r*7  [w^1)  +  v(1) /2] 
2p  n,  1  2c  2p 


Equations  (37)  and  (38)  are  the  same  as  those  of  the  Thau- 
Pao  solutions.  Equations  (39)  and  (40)  are  the  same  as  the  solutions 
presented  in  equations  (24)  and  (25) .  The  reason  we  obtain  different  forms 
of  particular  solutions  is  due  to  the  fact  that  the  following  expressions 
exist  to  generate  these  particular  solutions: 

(43)  (v2  +  k2)  [r*Vw^}+  k2r2  w^o)]  =  0 

%  %  ip  l 

i  =  1,  2. 

However,  although  we  can  find  the  different  forms  for  the  particular 

solutions,  after  satisfying  the  prescribed  boundary  conditions  for  some 

(2)  (2)  (2) 

specific  problems,  the  complete  solutions  w^  =  w^c  +  (i  =  1,2.) 

would  be  the  same  no  matter  which  forms  of  equations  (37)  to  (42)  are 


10 


chosen.  Similarly  different  expression  of  the  particular  solutions 
for  higher  orders  (n>2)  also  can  be  found.  Therefore,  here  we  shall 
designate  these  solutions  as  "Pseudo-particular  solutions"  instead 
of  "particular  solutions". 

As  for  the  perturbation  parameter  e,  equation  (9)  can  be 
rewritten  as 

[1  -  (v./v,)2] 

(44)  £  =  - 5 - 

2 [ 1  +  (v1/v2)^] 

In  an  elastic  solid,  as  Thau  and  Pao  pointed  out,  the 
numerical  value  of  the  parameter  £  is  limited  between  1/6  and  1/2. 

If  the  problem  of  SH-wave  in  two  different  welded  media,  or  acoustic 
waves  in  two  ideal  fluids  in  contact,  is  considered,  the  numerical 
value  of  e  is  1/2,  for  the  extreme  case  when  v^ — >  0  or  v^ — 

However,  there  is  no  limit  on  the  lower  bound  numerically.  The 
smaller  is  the  contrast  between  v^  and  v 2<  the  smaller  is  the  numerical 
value  of  e.  For  example,  for  the  case  of  two  shear  wave  velocities 
v^=3 . 7  km/sec,  v2=4.5  km/sec,  e  is  0.0966.  The  first  few  order  solu¬ 
tions  are  good  enough  for  a  small  quantity  of  e.  Therefore,  obtaining 
the  general  recurrence  formula,  no  matter  how  tedious  the  process  may 
be,  these  higher  order  perturbation  solutions  promise  to  solve  some 
long  standing  elastodynamic  problems  which  otherwise  are  analytically 


untractable. 


APPLICATION  OF  THE  SMS-VALUE  PERTURBATION  METHOD  TO 
SH-WAVES  IN  A  SLOPING  INTERFACED  HALF-SPACE 


The  geometry  of  the  problem  in  question  is  an  elastic 
wedge  overlying  an  elastic  bottom  as  shown  in  Figure  1.  Both  media 
are  assumed  to  be  homogeneous  and  istropic  solid.  The  inclined 
interface  is  assumed  to  be  perfectly  welded.  Restricting  a  line 
source  of  the  SH-wave  type  located  in  the  wedge,  we  can  consider  the 
problem  to  be  two-dimensional  so  that  the  particle  motion  of  SH- 
waves  is  confined  to  the  direction  of  z  only,  and  the  wave  propaga¬ 
tion  direction  is  in  the  (r,0)  of  the  polar  coordinate  system. 

The  overlying  wedge  medium  (  y ^ , p  ^ )  is  bounded  by  9=0  and  9=a;  the 
under-lain  bottom  medium  (^2' p boundec*  ^  6=ct  and  where 

and  p^,  i=l,2  are  the  rigidities  and  densities  of  the  overlying 
and  bottom  media  respectively.  A  line  source  of  SH-wave  is  located 

at  the  (r  ,0  )  in  the  overlying  wedge.  The  displacements  W. (r,9,t), 
o  o  1 

i=1.2.  satisfy  the  wave  equations  in  the  overlying  medium  and  bottom 
medium,  respectively 


(45) 


(46) 


2  ir 

r 


6 (r-rQ) 6 (6-60) 5 (t-0+) 


0 


where 


(47) 


6i  ”  ^j/Pi5 


i  =  1,2,  Bx<  32  . 
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The  conditions  to  be  satisfied  in  the  traction  free  boundaries  are 
sufficiently  described  by  the  vanishing  of  the  normal  components 
of  the  stress.  At  the  perfectly  welded  interface  between  two  dis¬ 
similar  media,  there  must  be  continuity  of  the  displacements  and 
stresses. 

By  taking  Laplace  transforms  of  the  equations  (45)  and 
(46)  with  respect  to  t,  the  new  transformed  equations  are 


2  TT 

(43)  (V2  -  h2)  w  (r,6,p)  - - 6  (r-r  )  6  (© -9  ) 

x  1  r  o  o 


(49)  (V2  -  h2)  w2(r,e,p)  =  0 


where 


(50) 


w  (r,0,p)  =  [W.  (r,9,t)]  =  /“  W.  (r,0,t)  exp(-pt)  dt 

1  -L  o  i 


and  hi  =  p/8i  i  =  1,2 


Using  the  present  RMS  perturbation  method,  we  obtained 
a  third-order  formal  perturbation  solutions  as 


(51) 


W.  =  W(0)  ♦  EW(1>  +  £2  w(2)  +  t3  W.(3) 

111  1  1 


where 


(52) 


2  2 

hI-h2 


2  (h2  +  h2  ) 


6l  -  el 


2  (fij  +  &22  ) 
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The  Zero-order  displacement  components  are: 


(53)  w 


(o) 


2  i 


n  -  *> 


ch  s(tt-0>+8<)  - 

-  sh  s(r>-0  )  sh  s9  + 
o 

+  ch  s9 

s 


K.  (hr  )  K .  (hr)  ds 
is  o  is 


(54)  w 


(o) 


2  i 


,  (o) 


2  '  _oo  -s  ch  S(1T  -0)  Kis(hro)Kis(hr)  ds 


with 


(55) 


-  y2  sh  s(-rr-a)[ch  s(ir-a+0o)  -  sh  s(ir-0o)  sh  sa  ]  + 
+  y  ch(ir-a)  [sh  s(ir-a+0o)  +  sh  s(it-0o)  ch  sa  ] 


[  y2  sh  s(ff-a)  ch  sa  +  y1  ch  s(u-a)  sh  sa  ] 


(56) 


,  (o) 


ch  sa  [sh  s  (ir-a+9  )  +  sh  s(tt-9  )  ch  sa  ]  + 
o  y 

+  sh  sa  [ch  s(T-a+9o)  -  sh  s(ir-9o)  sh  sa  ] 


[  y  sh  s(ir-a)  ch  sa  +  y^  Ch  s(ir-a)  sh  Sa  ] 


The  First-order  displacement  components  are: 


(57) 


(1) 


2  i 


ch  s(tt-9>+0<)  ” 

|-  sh  s(tt-90)  sh  s0  + 

|+  A(o)  ch  s9 
s 


**  hro  ^Kis+l(hro)  +  Kis-l(hr)]  Kis(kr) 


-  K  (hr)  r 
is  o 


3  Kia(hr) 

3  r 


ds 
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(58) 


i'¬ 


ll)  ^  ^  (1) 

w'  '  =  -  /  E  '  ch  s  ( tt— 6 )  K.  (hr)  ds  + 

2  _  2  -c0  S  IS 


IT 

2  i 


„  ,  ,  3K.  (hr) 

/  E  °  ch  s(ir-8)  K.  (hr  )  [-  r  vr^ -  ] 


ir 


2  -°°  S 


IS  o 


9  r 


with 


(59) 


E(1)  =  >5  hr  [K.  (hr  )  +  K.  .(hr  )]  E 

s  o  is+1  o  is-1  o 


(o) 

s 


The  Second-order  displacement  components  are  : 


(60)  w 


(2) 


2  i 


ch  s(7r-0>+6<)  - 


-  sh  s(ir-0  )  sh  s© 
o 

+  ch  s9 

s 


(61) 


(62) 


with 


T~  r  ro[Kis-l(hro)Kis+l(hr)+Kis+l(hro)Kis-l(hr)1  + 


+  (hr/2)  K  (hr  )[  K  (hr)  +  K  (hr)  ]  + 
is  o  is+2  is-2 

,2  2 

h  r 

+  - -  (  K.  ^  (hr  )  +  K.  .(hr  )  ]  K.  (hr) 

,  is+2  o  is-2  o  is 

2*2! 


(2)  (2)  A  (2) 

W-  —  w*  +  w 

2  2c  2p 


,<2) 

2c 


2  i 


; 

*00 


E^  ch  s(tt-6)  K.  (hr)  ds 
s  IS 


l 


ds 
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FUTURE  WORK 


The  formulation  of  the  higher  order  perturbation  method,  as 
an  example  for  solving  one  of  the  long  standing  wave  propagation 
problems,  i.e.  the  wedge  problem,  is  completed.  The  follow-up 
will  be  a  numerical  evaluation  of  the  formal  solutions  to  the 
problem  of  the  SH-wave  propagation  ir.  a  sloping  interfaced  nalf- 
ce 
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III. 


FINITE  ELEMENT  SOLUTION  TO  TRANSIENT  PROBLEMS 


OF  A  FINITE  CYLINDRICAL  CAVITY  SOURCE 


ABSTRACT 


This  paper  presents  the  numerical  solutions  to  the 
problem  of  the  transient  seismic  response  of  an  elastic  whole- 
space,  an  elastic  half-space,  and  a  layered  elastic  half-space 
due  to  a  finite  cylindrical  cavity  source  by  means  of  the  finite 
element  method,  based  on  the  principle  of  Virtual  Work. 

Synthetic  seismograms  of  both  the  radial  and  vertical 
displacements  of  the  transient  response  of  Model  I  -  a  whole- 
space,  Model  II  -  a  half-space,  and  Model  III  -  a  layer  half¬ 
space  for  a  finite  cylindrical  cavity  source  embedded  in  the 
medium  at  several  important  observation  points  are  obtained. 

Some  of  interesting  features  of  the  transient  response  are  as 
follows:  (1)  for  the  observation  points  inside  the  medium,  the 
reflected  waves  from  the  surface  boundary  of  the  half-space 
does  not  seem  to  have  significant  effect  on  the  amplitude  of 
the  wave  for  the  velocity  of  the  medium  assumed.  For  vertical 
displacements,  the  arrival  of  the  reflected  wave  is  identified 
as  a  step  in  time.  For  radial  displacements,  the  arrival  of 
the  reflected  wave  is  not  easily  identified,  (2)  the  amplitudes 
of  both  vertical  and  radial  displacements  at  the  observation 
point  on  the  surface  of  the  half-space  are  much  larger  than 
these  at  the  corresponding  observation  point  in  the  medium  of 
the  whole-space,  (3)  the  wave  forms  of  the  transient  response 
strongly  depend  on  the  dimensions  of  the  finite  cylindrical 
cavity,  i.e.  the  radius  and  length  of  the  cavity,  (4)  the  P  wave 
traveling  periodically  from  the  top  rim  to  the  bottom  rim  of 
the  cavity  along  the  vertical  interface  of  the  cylindrical 
'cavity  and  the  P  wave  periodically  oscillating  in  the  bottom 
(or  the  top)  face  of  the  cavity  with  a  motion  of  the  Bessel 
function  order  zero  type,  and  traveling  through  the  vertical 
face  of  the  cavity,  to  the  top  (or  the  bottom)  face  of  the 
cylindrical  cavity  also  with  a  motion  of  the  Bessel  function 
order  zero  type  are  dominant  in  all  three  models,  (5)  there  is 
no  vertical  displacement  at  any  point  along  the  middle  plane, 
whcih  are  symmetrical  with  respect  to  the  two  end  faces  of  the 
cylindrical  cavity,  in  Model  I  -  a  whole-space,  (6)  for  Model 
III  -  a  layer  half-space,  the  long  period  wave  is  predominant, 
traveling  periodically  along  the  verical  face  of  the  drilling 
hole  cavity  between  the  rim  of  the  contact  at  the  layer-interface 
and  the  top  rim  of  the  cylindrical  cavity,  (7)  the  weathering 
(low  velocity)  layer  acts  as  a  low  pass  filter,  all  the  high 
frequence  arrivals  are  filtered  out  in  the  seismograms  at  the 
observation  points  located  in  the  weathering  layer,  (8j  the  shape 
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'  of  the  input  forcing  function  is  we 11  preserved  on  each 

k  I  seismogram  of  the  vertical  displacement  but  not  on  the  radial 

displacement,  (9)  the  beat  phenomenon  is  observed  on  the 
radial  displacement  of  these  observation  points  located  in  the 
upper  region  of  the  cavity  source,  (10)  the  amplitudes  of  both 
the  radial  and  vertical  displacements  are  extremely  large  at 
those  observation  points  close  to  the  cavity  in  the  filled 
t  drilling  hole,  but  are  attenuated  rapidly  as  the  observation 

moves  away  from  the  cavity,  and  (11)  shear  waves  can  be 
generated  from  a  finite  length  cavity  source  with  only  normal 
stresses  applied  on  the  surfaces  of  the  cavity.  The  arrival 
of  the  shear  wave  are  clearly  identifiable  at  these  observation 
points  in  the  region  along  the  diagonal  line  of  the  cylindrical 
cavity. 
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1  INTRODUCTION 

For  the  transient  problem  of  a  layered  medium  due 
to  a  finite  cylindrical  source,  we  must  resort  to  numerical 
solutions,  such  as  the  finite  element  method,  which  provides 
a  mean  of  solving  a  class  of  transient  problems  of  elastic 
wave  scattering  and  diffraction.  For  this  reason,  we  have 
developed  the  finite  element  method  to  solve  transient 
problems  of  elastic  wave  propagation  in  a  layered  medium  due 
to  a  finite  cylindrical  source.  The  formulation  for  the 
finite  element  method  based  on  the  principle  of  virtual  work 
is  briefly  outlined  in  sections  2  through  4 .  It  should 
be  pointed  out  that  all  the  inherent  logarithmic  singularities 
are  being  removed  analytically  in  the  formulation.  Therefore, 
the  errors  due  to  the  singularities  will  not  be  introduced 
into  the  numerical  calculation.  The  applicability  of  the 
method  to  solve  practical  seismic  problems  has  also  been 
well  demonstrated  through  solving  a  series  of  the  transient 
problems  in  an  elastic  half-space  and  a  layered  half-space 
due  to  a  finite  cylindrical  cavity  source  as  shown  in  section 
5.  The  numerical  results  and  interpretation  are  presented 
in  section  6.  Conclusions  are  presented  in  section  7. 


2 


2.  DEVELOPMENT  OF  FINITE  ELEMENT  DISCRETIZATION  BASED  ON 
VIRTUAL  DISPLACEMENT  FOR  SOLVING  TRANSIENT  PROBLEM  OF 
A  FINITE  CYLINDRICAL  CAVITY  COMPRESSIONAL  SOURCE 


The  Principle  of  Virtual  Work  states  that  for  the 
equilibrium  of  an  elastic,  isotropic  and  homogeneous  body, 
the  total  internal  virtual  work  is  equal  to  the  total 
external  virtual  work  for  any  compatible  and  small  virtual 
displacements  imposed  onto  that  body  (Fung,  1965)  • 

f  x6et  dv  =  /  fb6Ufcdv  +  /  fs6Utds  1) 

V  V  s 


For  an  axial-symmetrical  geometrical  configuration 
with  a  dynamic  source  excitation,  we  adopt  the  cylindrical 
coordinates  system  (r,0,z).  The  quantities  in  1)  assume 
the  following: 


(a)  The  azimuthal-independent  displacement  field  Ub  and 
the  corresponding  strain  eb  are 

Ufc  =  [  ur,  0,  uzl 


=  ( 


rr 


'06' 


ezz ' 


0,  0, 


erzJ 


2) 


(b)  The  strain-displacement  and  the  corresponding 
stress,  and  the  stress-strain  relations  are 

£ij  ‘  <  ui,j  +  uj.i>/2 


=  [  T  ,  T „  -  ,  T  ,  0,  0,  T  ] 
rr  oo  zz  rz 


Tij  *  1  i£kkJij  +  2ucij! 


3) 


(c)  fS  and  fb 
respectively. 


are  the  surface  traction  and  body  force. 


1)  may  be  converted  to  a  system  of  ordinary 

differential  equations  by  the  following  arguments.  As 

the  volume  V  may  be  divided  into  n  subvolumes,  or  finite 

elements  Vm,  and  the  surface  S,  similarly  divided  into 

n  subsurfaces  S  ,  m  =  1,2,3, . . .  .N  We  then  have 

m 


N-*-  ® 


In  each  subvolume  (or  subregion)  V  ,  instead  of 

m 

determining  the  displacement  field  everywhere  in  Vm,  we 
find  the  displacement  vectors  at  only  a  discrete  finite 
number  of  points  in  V  ,  viz;  nodal  points,  under  the 
action  of  the  external  forces .  The  displacement  vector 
at  the  points  other  than  these  nodal  points  may  be 
determined  from  the  displacement  vectors  at  the  nodal  points 
through  some  interpolation  functions.  Within  each  element, 
the  displacement  field  um  may  be  expressed  in  the  form 


U111  =  NnS'Itl 


and  the  strain-displacement  relation  is 


ss 


where 

N111  *  a  known  matrix  of  the  spatial  interpolation 
function  for  the  mth  element, 

=  the  nodal  displacement  vector  of  the  mth  element, 
Bm  *  the  matrix  of  the  strain-displacement  relations 
of  the  mfc^  element. 


4. 


After  applying  equations  4) ,  5)  and  6) , 

equation  1)  becomes 


I  (M"Vm  +  K1”  +  Tm)  =0 
m=l 


7) 


with 


m"1  =  ! v m  dvm  8) 

k"*  =  /vm  (Bm)  dvm  9) 

Tm  -  /sm  (Nm)t  (fs)t  dsm  1Q 

where  M111 ,  Km,  Dm,  Tm  are  the  mass  matrix,  stiffness 
matrix,  matrix  of  constitutive  relation,  and  vector  matrix 
of  tractions  applying  on  the  surface  S,  respectively. 
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5. 


3.  INTEGRATION  OF  THE  INTERPOLATION  FUNCTION  ON  AN 
AXIAL-SYMMETRICAL  TRIANGULAR  RING  ELEMENT 

For  the  present  transient  problem,  which  possesses 
an  axial-symmetrical  geometry,  we  may  adopt  'solid 
toroidal  elements'  for  the  spatial  domain.  Since  the 
whole  formulation  is  independent  of  0,  the  integrals 
for  the  volume  and  surface  integrations  in  the  equations 
of  virtual  work  reduce  to  these  of  plane  surface  integrals 
and  line  integrals . 

Since  the  displacement  field  is  linear,  we  position 
the  element  arbitrarily  in  the  r-z  plane.  Constructing 
the  element  stiffness  matrix,  even  for  the  simplest 
case  of  a  linear  displacement  field,  we  must  take  into 
account  the  dependence  of  the  strain-displacement  equations 
on  the  inverse  of  the  radius,  particularly,  the  singularities 
along  the  axis  of  symmetry.  Applying  the  proper  interpola¬ 
tion  function  for  space,  and  using  the  pertinent  strain- 
displacement  relationship,  we  obtain  the  kernel  stiffness 
matrix  as  follows : 


in  which 


I1  =  IS  rdrdz 


I4  =//  drdz/r 


I 2  -IS  drdz 


I5  =//  zdrdz/r 


I3  -IS  zdrdz 


Ig  =//  z^drdz/r 


The  spatial  integrations  1^  through  Ig  have  carried 
out  analytically  so  that  the  inherent  logarithmic 
singularities  are  being  removed.  Otherwise,  the  errors 
due  to  these  singularities  will  be  introduced  into  the 
numerical  calculations.  The  terms  1^,  I2»  and  I3  yield 


I3  =  (1/6)  (r1+r2+r3)  [r-j^  (z2-z3)+r2  (z3-z1)+r3  (z1-z2)  ] 

l 

12  »  (1/2)  [r1(z2-z3)+r2(z3-z]L)+r3(z1-z2)  ] 

13  *»  (1/6)  (z1+z2+z3)  [r1(z2-z3)+r2(z3-z1)+r3(z1-z2)  ] 


The  terms  1^,  Ig,  and  Ig  involve  the  variable  r  in  the 
denominator  of  the  integrand  and  result  in  considerably 
complicated  expressions.  * 


8. 


Special  problems  arise,  when  the  joints  are  located 

on  the  axis  of  symmetry,  as  the  terms  with  ln(r^Aj)  and 

(r.-r.)  in  the  denominator  become  infinite.  In  this  case, 

x  j 

we  evaluate  the  integrals  by  L' Hospital's  rule. 


(ii)  In  the  case  r./r./r.  /r. ,  and  r.=0 

X  J  Jv  X  X 


I  =  [—2 — iS — i - iS — i — 1_]  in— X 

r.-r.  r. 

3  k  k 


I5  =  Y_.k  +  (1/4)  [(z^Zj)  {3zi+Zj)  +  (zk-2i)  (3zi+zk)] 

o  r  ■ 

-  (1/2)2.  In -I 

rk 


J6  *  z]k  +  “TT^  (ll2i2+5zi V22k2>  + 


2  .  *■  2  .  X*  • 

i-  j  (llz.2+5z.z.+2z.2)  -  (1/3)  z. 3  ln-JL 
18  1133  rk 


(iii)  In  the  case  ri=r k,  and  r^/r^O 


I4  -  t( - X.)  In— 1]  (zk-Zi) 

rj"rk  rk 


Ic  = - T\r  ■  z, +r  .  z.-2r.r.  z.)  In — i- 

5  2(ri-rk)3  3  k  ]  i  3k:  ^ 


zk'zi 


— ^[22  .  (rk+r .)-  Orj-r^)  ( Zfc*z  . )  1 


•j  ~k 


^Zk~2i^ri  22  2 

I*  = - — 3_  [r j~  (z.  +z.  z, +z  .  )  +3r.  z  .  (r,  z  .-r  .z. - 

6  3(r.-rk)3  3  k  i  k  i  k  j  k  3  3  k 


^Zk~ Zi^  2  2 

r.z  )]  In  +  - — — — y  [z.2(3r.  -15r  .r  - 

]  1  rk  18(rj-rk)^  3  3  3  k 


6rk2)+zj(2k+zi><6rj2+15rjrk-3rk2)  -  <zk2+zi2+ 


zizkHUrj2-7rjV2rk2)) 


(iv)  In  the  case  r.=r.^r.  ,  r.=r.^0,  and  r.=0 

*  ^  J  IK  j 


I4  =  -(zk-zi) 


V*i 

!5 - (2zj+z.+zk) 


z  ^  z 

I  - k — i>  (6z.2+2z,_2+2zJ2+3z.zJ  +  3z  .z, +2z.z,  ) 


(v)  In  the  case  when  r.=r.=0,  and  r,=0,  it  can  be  shown 

1  j  JC 

that  the  terms  containing  1^,  1^ ,  and  1^  do  not  appear 

in  the  stiffness  matrix  that  the  displacements  u.=u.=0. 

i  3 

All  these  integrals  can  be  used  as  a  basis  for  deriving 
the  similar  integrals  of  nonlinear  interpolation  functions, 


4.  TIME  INTEGRATION  SCHEME 


After  careful  and  critical  tests  and  evaluation  of 
various  three-dimensional  finite  element  formulations, 
it  was  found  that  the  combination  of  lumped  mass  matrix 
representation  and  explicit  central-differences  time 
integration  of  three-dimensional  finite  element 
formulation  appears  to  be  the  best  in  solving  transient 
problems  of  elastic  wave  propagation.  An  important 
advantage  of  using  a  lumped  mass  matrix  is  that  the  matrix 
is  diagonal,  and  as  it  will  be  seen  later,  the  numerical 
operations  for  obtaining  the  solution  to  the  dynamic 
equation  of  equilibrium  are  simplified. 

The  scheme  of  time  integration  for  the  dynamic 
structural  problem  with  the  equation  of  motion 


Mtf  +  KU  =  F 


17) 


can  be  considered  as  a  subclass  of  the  explicit  methods  of 
the  form 


11. 


N 

U  . ,  +  Z  (a.U  .  +  b .KU  . )  =  0 
n+1  j_Q  3  n  3  3  n-3 


with  N=2,  aQ=-2 ,  a^=l,  bg=At  ,  b^=0. 


Therefore,  we  obtain 


18) 


U(t+At)  =  U(t)  +  0  (t) At 

C (t+At)  =  0(t)  +  M_1KU (t+At) At  19) 


which  is  the  central  differences  scheme.  Krieg(1970) 
and  Baylor  et  al(1974)  have  proved  that  the  central 
differences  scheme  has  the  largest  stability  region. 
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5.  TRANSIENT  PROBLEM  DUE  TO  A  FINITE  LENGTH  CYLINDRICAL 
CAVITY  IN  AN  ELASTIC  HALF-SPACE  AND  LAYERED  HALF- SPACE 


In  this  dissertation,  with  a  view  to  simulate 
conditions  in  the  neighborhood  of  a  drilling  hole,  the 
source  chosen  is  a  cylindrical  cavity  of  length  %  =  6 
feet,  and  radius  a  =  0.5  feet,  centered  at  the  origin  of 
coordinate  and  embedded  in  an  infinite  medium.  Figure  1. 

For  comparison,  we  simulate  the  second  model  with  the 
half-space,  the  surface  of  the  half-space  is  located  at 
25  feet  above  the  top  end  of  the  cylindrical  cavity,  as 
shown  in  Figure  2.  The  half-space  medium  is  character¬ 
ized  by  the  elastic  constants  E,  a,  p  as  medium  2  given 
in  Table  1. 

In  seismic  prospecting,  a  layer  occurs  at  the  surface 

of  the  earth  which  is  unconsolidated,  and  often  heterogeneous 

with  low  wave  velocity,  so  called  "weathering  layer". 

The  thickness  of  the  weathering  layer  may  vary  from  almost 

zero  to  several  hundred  feet;  however,  the  most  common 

thickness  is  from  about  25  feet  to  100  feet.  For  our 

computation,  we  simulate  the  weathering  layer  with  a 

thickness  of  20  feet.  The  cylindrical  cavity  is  embedded 

5  feet  below  the  weathering  layer,  including  the  case  of 

a  filled  drilling  hole.  The  geometrical  setting  is  shown 

in  Figure  3.  The  medium  of  the  weathering  layer  and 

the  filled  drilling  hole  are  characterized  by  the  elastic 

constants  E^ ,  and  E^ ,  respectively,  as 

given  in  Table  1.  The  finite  element  meshes  are  formed 

according  to  the  dimensions  as  described  above  and  shown 

* 

in  Figure  4^5  Model  I  has  687  nodal  points  and  1268 


25  feet 


20  feet 


5.  The  Geological  Setting  for  Model.  III. 


Figure 


elements  and  Model  II  and  III  have  729  nodal  points 
and  1348  elements.  The  size  of  the  finite  element 
models  is  chosen  in  such  a  way  that  the  artificial 
reflections  will  not  reach  the  observation  points  within 
the  time  interval  of  0  -  15  milliseconds. 

The  system  is  set  into  motion  by  the  application  of  a 
sudden,  normal  stress  applied  uniformly  on  both  the 
vertical  walls  and  the  two  ends  of  the  cylindrical  cavity. 

The  time  variation  of  the  applied  stresses  are  assumed  to 
be  given  by  a  modified  Heaviside  function  as  shown  in 
Figure  6.  The  rise  time  of  the  pulse  is  assumed  to  be 
16At,  where  At=20  ^seconds,  is  in  order  to  satisfy  the 
Courant's  stability  condition. 

In  order  to  study  the  transient  response  in  the  region 
along  the  vertical  surface  of  the  cylindrical  cavity,  the 
observation  points  along  LINE  A,  B,  C  are  selected  (see 
Figure  7) .  In  order  to  study  the  diffraction  behavior 
of  the  rims  of  the  cylindrical  cavity  and  the  comparison 
between  the  transient  response  in  the  weathering  layer  and 
the  half-space,  we  choose  the  observation  points  to  be  located 
along  LINE  D  and  E.  The  observation  points  along  LINE  F, 

G, H, I , J, K  are  chosen  to  investigate  the  transient  response 
in  the  upper  region  of  the  cylindrical  cavity  and  in  the 
filled  drilling  hole.  For  the  surface  response  as  encon- 
tered  in  seismic  exploration,  the  observation  points  are 
selected  to  be  located  along  LINE  L. 

For  the  limited  computer  time  available  to  the  present 

study,  only  the  computational  results  for  the  displacement 

u  and  U  in  the  time  interval  of  0  -  15  milliseconds  for 
r  z 

the  observation  points  along  LINE  A  through  L  of  Model  I, 

II,  III  are  obtained  as  presented  in  Figure  10  through 

21. 
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Figure  7a.  The  Locations  of  the  Observation  Points 
in  the  Numerical  Solution  for  Model  I. 
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re  7b.  The  Locations  of  the  Observation-Points 
ie  Numerical  Solution  for  Model  II.  * 


LINE  F,C,H  LINE  0 

L 

r 

Xi'Ul'°l 

r 

L 

/ 

I  A72 

214  215  216  217  218  219  220 

253  254  255  256  257  258  259 

293  294  295  296  297  298  299 

V 

r 

r 

\ 

•400 

V 

r 

LIME  I.J.K  trN*  E 


X2'VE2 


7c.  The  Locations  of  the  Observation  Points 
Numerical  Solution  for  Model  III. 


6.  INTERPRETATION  THE  NUMERICAL  RESULTS 


A.  Ray  Paths  for  Model  III 

It  appears  to  be  instructive  to  use  the  simple 
concept  of  ray  to  show  the  various  arrivals  of  waves 
in  Model  I,  II,  and  III.  As  the  ray  paths  for  Model 
I  and  II  are  relatively  easily  traceable,  in  order  to 
avoid  duplication,  we  choose  to  show  schematically  the 
principal  ray  paths  for  Model  III  of  a  finite  cylindrical 
cavity  source  located  in  a  two-layered  half-space  with 
a  filled  drilling  hole  with  unconsolidated  material. 

It  is  assumed  that  the  upper  layer  models  a  weathering 
layer,  whcih  is  underlain  by  a  relatively  consolidated 
half-space.  The  filled  drilling  hole  is  properly  tapped, 
and  there  are  no  blowouts.  The  cavity  simulates  a 
cylindrical  charge.  The  following  notations  for  the 
expected  arrivals  of  waves  are  adopted: 

(1)  Subscripts  1,  2  and  3  refer  to  wave  velocities 
in  the  upper  layer;  i.e.  the  weathering  layer,  the  under¬ 
lying  half-space  and  the  filled  drilling  hole,  respectively 

(2)  Superscript  t  refers  to  the  rim  of  the  filled 
drilling  hole  on  the  surface,  the  rim  of  the  contact  at 
the  layer- interface,  the  top  rim  of  the  cylindrical  cavity, 
and  b  refers  to  the  rim  of  the  contact  at  the  layer 
interface,  the  bottom  rim  of  the  filled  drilling  hole,  and 
the  bottom  rim  of  the  cylindrical  cavity,  so  that  for 
example,  lb  =  2t,  2b  =  3t. 

(3)  Subscripts  J  refers  to  cylindrical  waves  on  the 
top  and  bottom  faces  of  the  filled  drilling  hold  and  the 
cylindrical  cavity  of  the  Bessel  function  order  zero  type. 


i 


t 


{. 


24. 

(4)  Superscripts,  IN,  2N,  and  3N,  refer  to  the  repeated 
multi-ray  paths  along  the  surface  of  the  filled  drilling 

hole  between  It,  lb;  2t,  2b;  and  along  the  cavity  respectively. 
N  is  the  number  of  repeated  multi-ray  path. 

(5)  Primes,  refer  to  refracted  waves. 

The  ray  paths  of  the  principal  wave  arrivals  are 
shown  in  Figure  8  and  Figure  9 ,  for  P-waves ,  and  for 
P  converted  into  S,  respectively.  The  ray  paths  for  S 
waves  and  for  S  converted  into  P  waves  can  be  similarly 
traced,  according  to  the  Snell's  law.  The  following  are 
the  identification  of  various  waves  for  Models: 

(1)  Direct  Waves, 

p,  S:  P  and  S  waves  directly  from  a  point  of  the  cavity 
source  to  the  observation  point. 


(2)  Direct  Diffracted  Waves, 

P3t,  P3b,  S3t,  S3b:  Diffracted  P  and  S  waves 
traveling  from  the  top  and  bottom  rim  of  the  cylindrical 
cavity  source. 


P3bP 
*2  P2' 


P3bP 
P2  P2 


P3tS 
2  S2 ' 


p3bq  . 
P2  S2* 


(3)  Diffracted  and  Reflected  Waves, 

Diffracted  P  wave 
traveling  from  the  top  and  bottom  rim  of  the  cylindrical 
cavity  source  to  the  layer  interface  and  reflected  as  P 
wave  or  converted  as  S  wave,  respectively. 

Diffracted  S  wave 
traveling  from  the  top  and  bottom  rim  of  the  cavity  source 
to  the  layer  interface  and  reflected  as  P  wave  or  converted 
as  S  wave,  respectively. 


S3bP 
S2  P2 ' 


S3bp 
S2  P2 ' 


S3tS 
S2  S2 ' 


S3bS  • 

s2  s2. 
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(4)  Refracted  Waves, 


p2  pi'  p2  pi#  p2  si'  p2  si:  p  wave  travelin9 

from  the  top  and  the  bottom  rims  of  the  cavity  source  to 
the  layer  interface  and  refracted  into  the  first  medium  as 
P  wave,  or  converted  as  S  wave,  respectively. 


S2  Sl'  S2  Sl'  S2  pi'  S2  Pi:  S  wave  travelin9 

from  the  top  and  the  bottom  rim  of  the  cavity  source  to 
the  layer  interface  and  refracted  into  the  first  medium  as 
P  wave  or  converted  as  S  wave,  respectively. 


3t  3 1 

P^,  P2  S3:  P  wave  travelin9  vertically  from  the 
top  rim  of  the  cavity  source  along  the  interface  between 
the  filled  drilling  hole  and  the  layered  medium,  with  the 
velocity  of  the  half-space  and  refracted  into  the  filled 
drilling  hole  as  P  or  converted  to  S. 

2 1  2 1 

P^  P^,  P^  S^:  P  wave  traveling  vertically  upward 
from  the  rim  of  the  contact  between  the  layer  interface  and 
the  filled  drilling  hole  and  refracted  into  the  filled 
drilling  hole  as  P  or  converted  to  S. 

3t  3t 

S2  ,  Sj  P^:  S  wave  traveling  vertically  from  the 
top  rim  of  the  cavity  source  along  the  interface  between 
the  filled  drilling  hole  and  the  layered  medium,  with  the 
velocity  of  the  half-space  and  refracted  into  the  filled 
drilling  hole  as  P  or  converted  to  S. 

2t  2t 

S1  S3'  ®i  P3:  s  wave  traveling  vertically  upward 
from  the  rim  of  the  contact  between  the  layer  interface  and 
the  filled  drilling  hole  and  refracted  into  the  filled 
drilling  hole  as  P  or  converted  to  S. 


(5)  Multi-Diffracted  Waves, 
pfpf*  p23tpfpfpf<=  ?23IIIP23b> , 

pfpf,  p3bP3tp3bpf<=  P3IIIpfi, 


.  ^  jaht  ■  -f  a  u  itj-fli. 
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P 

P 


3tP3bs 
2  P2  S 


3bP3ts 
2  P2 


2' 

2' 


p3tp3bp3tp3bq 
2  2  2  2  2{ 

P3bp3tP3bP3tS2(= 


P 

P 


3lIIc3b 
2  S2 

3IIIc3t 
2  b2 


)  , 


)  , 


Diffracted  P  wave  traveling  periodically  from  the  top 
rim  to  the  bottom  rim  of  the  cavity  along  the  vertical 
interface  of  the  cylindrical  cavity  and  diffracted  as 
P  wave  or  converted  as  S  wave  respectively  at  the 
top  (or  bottom)  rim  of  the  cavity. 


p2bp2t  p2bp2tp2bp2t  f_  p2IIIp2t, 

*2  *2  '  P2  P2  p2  *2  P2  P2  '  ' . 

2b  2t  2b  2b  2t  2b  2t  2b,  p2IV  2b. 

P2  2  P2  '  2  2  P2  P2  P2  P2  P2  '  ' . 

P2b  2t  2b  2t  2b  2t  2III  2t. 

2  2  s2  '  2  2  2  P2  *2  P2  &2  ;  . . 

2b  2t  2b  2b  2t  2b  2t  2b  P2IVS2t^ 

P2  P2  P2  S2r  P2  P2  P2  P2  P2  S2(_  P2  S2  } ' 


Diffracted  P  wave  traveling  vertically  along  the  hole 
between  the  top  of  the  cylindrical  cavity  and  the  layer 
interface  and  diffracted  as  P  wave  or  converted  as  S 
wave  respectively  at  the  top  rim  of  the  cylindrical  cavity 
or  the  layer  interface. 


(6)  Special  Waves, 


vMVf' 


PJP23kp2tpjP2t' 


pjp2tp2bpjp2bp2tpjp2t(=  PJP2IIp23t|  ’ 


P  wave  traveling  periodically  from  the  bottom  (or  the  top) 
face  of  the  cavity  with  a  motion  described  by  the  Bessel 
function  order  zero,  through  the  vertical  face  of  the  cavity, 
to  the  top  (or  the  bottom)  face  of  the  cylindrical  cavity 
also  with  a  motion  described  by  the  Bessel  function  order 


zero. 


B.  Detailed  Interpretation  of  Synthetic  Seismograms 

Figure  10  through  21  show  synthetic 
seismograms  of  both  the  radial  and  vertical  displacements 
U  and  U  of  the  transient  response  of  Model  I  -  a  whole 

i  Z 

space,  Model  II  -  a  half-space,  and  Model  III  -  a  layer 
half-space  due  to  a  finite  cylindrical  cavity  source 
embedded  in  the  medium  at  the  observation  points  along 
lines,  A  to  L,  as  shown  in  Figure  7. 

We  here  shall  attempt  to  identify  the  principal  wave 
arrivals,  by  means  of  both  the  recognition  of  the  change 
of  the  characteristics  of  the  wave  forms  and  the  expected 
arrival  times  of  these  principal  waves  based  on  the  ray 
tracing  technique. 


Line  A 


Figure  10  shows  synthetic  seismograms  of  Line  A, 
at  the  nodal  points  214,215,216,217,218,219  and  220  for 
Model  I,  II,  and  III.  (See  Figure  7  for  the  locations 

of  the  nodal  points)  In  Figure  10  a^,  the  radial 

displacement  U 


the  arrivals  of  each  trace,  P2,  S2b, 


PJP2bp2tpjP2'  P2tp2bp2t'  and  the  multiPles  of  PjP2bp2tpjP2' 
P2tp2bp2b  are  identified.  The  paths  of  these  arrivals 
are  shown  in  Figure  8  and  Figure  9.  In  addition  to 
these  arrivals,  by  comparing  the  radial  and  vertical 
displacements  of  Model  I  with  these  of  Model  II,  the  arrival 
of  the  wave,  P^P^  identifie<3  as  a  step  in  time  and  is 


indicated  in  Figure  10  a2 

p  p  p^kp^p  p  _3tp3b  3t 
2'  *\r2  *2  PJP2 '  P2  p2  P2  ' 

P^P^S. 


1  •  The  arrivals  of 

p2bp2tp  p2bp2t_  p2b  2t_ 

V2  *2  P2  *2  P2  S2'  P3  P3  P2' 


pMv 

v2  F2  P2 


S°S  tp,  and  the  multiples  of  P^P^P.,  P‘“p‘ws^, 
Sv,  T4.  4.  I  2  2  2 

•  S"d  PJP2  P2  PJP2 


,2b„2t„  .,2b  2 1, 

2  P2  * 

are  similarly  identified  in 
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Figure 


10  a3  for  Model  III. 


As  the  P  wave  velocity  in  the  half-space  is  assumed 
to  be  5000  ft/sec,  the  travel  time  of  1.2  milliseconds  for 
a  P  wave  is  equivalent  to  the  travel  distance  of  6  feet, 
which  is  the  length  of  the  cylindrical  cavity.  In  Figure 


10  an ,  the  periodic  wave  train  immediately  following  the 
arrivals  of  ^2  P2  ?2  with  a  period  of  about  1.2  milli¬ 
second  marked  by  "A"  on  the  seismograms  begining  at  4 . 2 


milliseconds  with  respect  to  the  origin  time,  appears  to 


fit  the  traveling  path  of  the  type  of  P  wave  P^P^p^t 


2  ‘  2  2 

traveling  periodically  from  the  top  rim  to  the  bottom  rim 
of  the  cavity  as  P  wave  along  the  vertical  interface  of 
the  cylindrical  cavity.  The  wave  train  with  a  period  of 


1.6  milliseconds  marked  by  "B"  on  the  trace  214  in  Figure 


10  a1  is  indicated  as  PjP^P^PjP^k,  as  the  travel  time 


of  1.6  milliseconds  for  P  in  the  half-space  is  equivalent  to 
the  traveling  distance  of  8  feet.  The  P  wave  PjP^P^PjP^13 
apparently  is  traveling  back  and  forth  periodically  around 
the  cylindrical  cavity.  Because  the  distances  between 
-each  of  the  observation  point  and  the  top  and  bottom  rims 
of  the  cylindrical  cavity  are  different,  waves  originated 
from  the  top  rim  and  from  the  bottom  rim  of  the  cylindrical 
cavity  interfere  constructively  and  destructively ,  depending 
on  the  distances  between  the  observation  point  and  the  top 


and  bottom  rims  of  the  cavity.  In  Figure 


.10  a_ ,  for  the 

velocity  of  the  medium  assumed,  reflected  waves,  P^  Pj^  from 


the  surface  boundary  of  the  half-space  does  not  seem  to  have 


significant  effects  on  the  amplitude  of  the  wave,  compared 
3t 
2 


with  P^v"  in  a  whole  space,  as  in  Figure  10  a,.  Seismo- 


ll* 


grams  in  Figure  10  a^,  and  especially 


10  a^ ' ,  show  the 


illffUfiltiSfliiriMiHi 
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superposition  of  the  short  period  waves  of  about  1.2  milli¬ 
seconds  on  the  long  period  wave  of  about  5  milliseconds,  as 
marked  by  "L".  The  long  period  waves  observed  on  these 
seismograms  are  essentially  generated  in  the  filled  drilling 
hole  with  a  very  low  velocity  as  radiated  cylindrical  waves, 
as  the  refracted  waves  from  the  weathering  layer  arrive 
later  than  the  time  of  1.5  milli-second.  The  travel  time 
of  the  P  wave  in  the  filled  drilling  hole  from  the  top  of 
the  cylindrical  cavity  to  the  layer  interface  is  5  milli¬ 
second,  corresponding  to  a  distance  of  5  feet.  The  travel 
distance  of  the  wave  pf°P^P_  is  exactly  5  feet.  Therefore, 
these  long  period  waves  are  identified  as  P^  P^  ,  traveling 
along  the  hole  vertically  between  the  top  of  the  cylindrical 
cavity  and  the  layer  interface.  The  input  step  function 
is  not  as  clearly  preserved  in  the  radial  displacement  in 
Figure  10  a^  a2 ,  and  a^  as  in  the  vertical  displacements  as 
shown  in  Figure  10  a|,a2  and  a^ •  For  a  finite  cylindrical 
cavity  source,  the  preservation  of  the  input  step  largely 
depends  on  the  configuration  of  the  observation  point. 

All  the  observation  points  in  Figure  10  a^,  a2,  and  a^ 
are  located  along  a  line  perpendicular  to  the  axis  and 
parallel  to  the  top  and  bottom  faces  of  the  cylindrical 
cavity.  The  normal  applied  stresses  on  the  vertical 
face  of  the  cavity  contribute  to  the  radial  displacement 
in  a  radially  outward  direction,  but  the  applied  stresses 
on  the  end  faces  of  the  cavity  tend  to  elongate  the  cavity 
in  an  axial  direction  and  contribute  to  the  radial 
displacement  in  a  radially  inward  direction  along  the 
vertical  face  of  the  cavity.  The  combination  of  the 
applied  stresses  on  the  both  end  faces  and  the  vertical 
face  of  the  cavity  gives  the  resultant  complicated  radial 
displacement.  Because  Nodal  Point  214  is  only' 2.5 

feet  from  the  top  rim  of  the  cylindrical  cavity,  the 
contribution  due  to  the  applied  stresses  on  the  both  end 
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faces  of  the  cylindrical  cavity  to  the  radial  displacement, 

at  Nodal  Point  214  is  approximately  the  same  order  of 

magnitude  as  that  due  to  the  applied  stresses  on  the  vertical 

face  of  the  cavity,  as  reflected  through  a  radially  outward 

3t 

response  of  the  wave  forms  between  the  arrivals  P0  and 

3b  ^ 

S2  •  As  the  distances  between  Nodal  Points  215  to  220, 

and  the  cylindrical  cavity  increase,  the  contribution  due  to 

the  applied  stresses  on  the  both  end  faces  of  the  cavity 

becomes  pre-dominant,  as  reflected  through  the  radially 

3 1 

inward  response  of  the  wave  forms  between  the  arrivals  P_ 
and  ,  as  designated  by  "C" .  In  Figure  10  a^ ' ,  and 
a2 *  t  the  vertical  displacements  U  ,  the  wave  arrivals  P2» 

S2b'  P2tp2bp2ts2'  PJP2IIIp2t'  P^P^P^P^hp^  and  the  multiPles 

of  p2tp2bp2ts2'  and  P2tp2bp2tp2t>P2t  are  identified* 

The  ray  paths  of  these  arrivals  are  given  in  Figures  8 

and  9.  In  addition  to  these  arrivals,  the  arrival  of  a 
3t 

step  P2  P2/  in  Model  II  Figure  10  a2 ' ,  is  identified. 

In  Figure  10  a  '  of  the  vertical  displacement  U  of  Model 
III,  the  arrivals,  P2,  P^P^,  P2bp2ts2'  P^P^*  P^P^^' 

and  P3bp3tS2  and  the  multiples  of  P^P^/  P^P^Sj,  P^P^ 

p2bp2tS2  are  identified.  The  wave  arrival  of  P^P]^^^ 

predominant  on  all  the  seismograms  in  Figure  10  a,’  , 

&2 1  and  a3  ’  •  The  arrival  of  P^  Pj  Sj  contrarily  is  only 

predominant  on  the  seismograms  of  Figure  10  a,’  .  The 
3b  3t  3t 

wave,  pjp2  P2  pjP2  '  identified  in  each  seismogram  of 
Figure  10  a^ '  and  a2 ' •  It  has  a  period  of  about  1.6 
milliseconds ,  as  specifically  marked  "A"  on  the  trace  of 
Nodal  Point  217  in  the  time  interval  of  8.3  milliseconds 


and  9.9  milliseconds, 


The  identification  of  the  wave  is 


based  on  the  fact  that  the  P  wave  velocity  in  the  half-space 
is  5000  ft/sec,  the  travel  time  of  1.6  milliseconds  for 
the  P  wave  is  equivalent  to  the  travel  distance  of  8  feet. 


] 
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Therefore,  the  wave  PjP2  PjP2  is  the  type  of  the  wave 
traveling  periodically  from  the  bottom  (or  the  top)  face 
of  the  cavity  with  a  motion  described  by  the  Bessel  function 
order  zero,  through  the  vertical  face  of  the  cavity,  to 
the  top  (or  the  bottom)  face  of  the  cylindrical  cavity  also 
with  a  motion  described  by  the  Bessel  function  order  zero. 
The  ray  path  of  this  type  of  wave  is  shown  in  the  following 
diagram. 


Line  B 


Figures  11,  give  synthetic  seismograms  of  the 
radial  displacement  Ur  and  vertical  displacement  Uz  of 
Line  B  at  Nodal  Point  253,254,255,256,257,258  and  259 
for  Model  I,  II,  and  III  as  shown  in  Figure  7.  The 
waves  P. 


p3t 


P3tP3b 
P2  P2  ' 


p3b  3t 
2  2  ' 


PJP2tp2bpjP2b  3nd  the 


multiples  of  P^P^' 
identified  in  Figure  11  a^.  Model  I,  the  whole  space. 


E’jp2tp2bpjp2b  are  oil 


The  reflected  and  converted  waves  from  the  layer  interface 
and  the  diffracted  from  the  top  rim  of  the  filled  drilling 

hole,  2t .  ^2^2 '  P2bp2'  P2ts2'  P2bs2'  ^^2^^ 1 

PJP2tp2bpjP2b  and  the  multiPles  of  p2tp2bp2fc'  P2tp2bp2tp2tf 
P^P^P^kp  are  identified  in  Figure  11  a,  for  Model 


Pj?2  P2  pjP2  are  i-et^  in  Figure  11  a^  for  Model 

III,  The  predominant  waves  observed  include:  (1)  P^P^P^' 
traveling  periodically  from  the  top  to  bottom  along  the 
vertical  face  of  the  cylindrical  cavity  with  a  velocity  of 
5000  ft/sec  and  a  period  of  about  1.2  milliseconds.  (2) 
PjP^P^kpj-^k  *  travels  periodically  from  the  bottom 
Cor  the  top  )  face  of  the  cavity  with  a  motion  described 
by  the  Bessel  function  order  zero,  through  the  vertical 
face  of  the  cavity,  to  the  top  (or  the  bottom)  face  of  the 
cylindrical  cavity  also  with  a  motion  described  by  the 
Bessel  function  order  zero. 


The  normal  applied  stress  on  the  vertical  face  of  the 
cavity  contributes  to  the  radial  displacement  in  a  radially 
outward  direction,  (i.e.  the  sense  of  direction  on  the 
seismogram  is  "positive")  but  the  applied  stresses  on  the 
end  faces  of  the  cavity  tend  to  elongate  the  cavity  in  an 
axial  direction  and  contribute  to  the  radial  displacement 
in  a  radially  inward  direction  (i.e.  the  sense  of  direction 
on  the  seismogram  is  "negative")  along  the  vertical  face 
of  the  cavity  everywhere.  Morever,  there  are  stress 
concentrations  around  both  the  top  and  the  bottom  rims  of 
the  cavity.  As  the  location  of  the  Nodal  Point  253  is  only 
2,5  feet  from  the  vertical  face  of  the  cavity  and  an  equal 
distance  of  3.9  feet  from  both  the  top  rim  and  the  bottom 
rim  of  the  cavity,  the  first  arrival  P  wave  from  a  minimum 
distance  between  the  observation  point  Nodal  Point  .253  and 
the  vertical  face  of  the  cavity  contributes  to  the  radial 
displacement  in  a  radially  outward  direction;  the  diffracted 
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arrivals  P^  and  P^b  from  both  the  end  faces  of  the  cavity 

contribute  to  the  radial  displacement  in  a  radially  inward 

direction.  Subsequently,  the  arrivals  of  the  diffracted 

waves  P^P^'  and  anc^  t^ie^r  multiples  contribute  to 

the  radial  displacement  to  have  an  oscillatory  appearence. 

Therefore  the  arrivals  of  the  waves,  P,  P?*",  P^,  P^P-h, 

3b  3t  ^  ^  c.  * c 

P2  P2  9-*-ve  an  impulsive  like  wavelet  in  the  time  interval 

of  0  and  3  milliseconds,  and  the  arrivals  of  the  multiples 

of  the  diffracted  waves  P^P^  and  P^P^^ivs  the  resultant 

oscillatory  radial  displacement  and  followed  by  a  long 

period  wave  train  as  clearly  shown  in  the  seismograms  of 

Nodal  Point  253  of  Model  I,  II,  and  III  (see  Figure  11 

al'  a2 '  and  a3* ‘ 

As  expected,  there  are  no  vertical  displacement  at 
every  nodal  point  along  Line  B,  symmetrically  with  respect 
to  the  two  end  faces  of  the  cylindrical  cavity  in  Model  I, 
in  the  whole  space  as  shown  in  Figure  11  a^ ' .  Because 
the  applied  stresses  on  the  both  end  faces  of  the  cavity 
are  equal  and  in  an  opposite  direction,  the  resultant 
vertical  displacement  along  Line  B  is  therefore  zero;  Like¬ 
wise  the  stresses  applied  to  the  vertical  face  of  the  cavity 
are  in  the  r-direction,  the  resultant  vertical  displacement 
along  Line  B  is  also  zero.  In  Figure  11  a2 1  ,the  vertical 
displacement  along  Line  B  for  Model  II,  a  half  space,  is 
zero  until  the  arrival  of  P2tp2'  ^*e'  p  wave  is  reflected  at 
the  surface  of  the  half-space.  The  waves  for  Model  III  are 
complex  as  shown  in  Figure  a^’ .  There  are  eight  principal 
arrivals  accounted  for,  e.g.  P^P^f  P2tp2'  P2ts2'  S2bs2' 

S2 ts2  *  PobP9tp7bp?ts->  *  P^P^P,  and  P^P^S?  ,  and  the 


2  r 2  3  3  02 


3  3  2 


3  3  s2 


multiples  of  p2bp2tp^  and  p2bp2ts^ 


Waves  with  a  period  close  to  5  milliseconds  marked  by 
"A"  are  clearly  identifiable  on  virtually  all  the  seismograms 
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in  Figure  11  a3'  ;  they  are  p3  P3  P?  traveling 
periodically  along  the  vertical  face  of  the  drilling  hole 
between  the  two  rims  at  2b  and  2t  with  the  velocity  of 
the  filled  drilling  hole. 


Line  C 


Figure  12  gives  synthetic  seismograms  of  the 
radial  displacement  U  and  the  vertical  displacement  Uz 
of  Line  C  at  Nodal  Point  293,294,295,296,297,298,  and 
299.  Line  C  lies  in  the  plane  of  the  bottom  face  of  the 
cavity,  as  Line  A  lies  in  the  plane  of  the  top  face  of  the 
cavity.  The  locat '  ns  of  these  nodal  points  are  shown  in 
Figure  7.  As  expected,  the  radial  displacements  of 
all  these  nodal  points  as  shown  in  Figure  12  a^  are 
identical  to  these  as  shown  in  Figure  10  the  vertical 
displacements  of  all  these  nodal  points  as  shown  in  Figure 
12  a^,  are  also  identical  to  these  of  Figure  10  , 

except  the  vertical  displacements  of  Line  C  are  in  an 
opposite  direction  of  these  of  Line  A.  The  seismograms 
of  the  radial  and  the  vertical  displacements  as  shown  in 

Figure  12  a„  and  a_  *  are  identical  to  these  as  shown  in 

*■  3t 

Figure  12  a^  and  a^',  until  the  arrival  of  P2  P2  anc^ 

p2bP2  as  in  the  case  of  Line  A.  In  comparison  of  the 

radial  and  the  vertical  displacements  Line  C  in  Figure 


11-12  a3  and  a3'  with  these  of  Figure  12  a.  and  a^' 
the  arrivals  of  P2tp2'  P2ts2'  P2^S2'  P3^P3tp2' 


P2bP2tS 
P3  P3  S2 ' 


P2bP2tP2bP2  and  P2bP2tP2bp2  are  easily  identified 


on  the  seismograms  of  Model  III  as  shown  in  Figure 


a3  and  a3  ' . 
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Line  D 


Figure  13  shows  synthetic  seismograms  of  the  radial 
displacement  Ur  and  vertical  displacement  Uz  of  Line  D  at 
Nodal  Point  193,172,152,131,111,90  and  70  for  Model  I,  II, 
and  III.  As  shown  in  Figure  7,  Line  D  is  situated 
diagonally  with  respect  to  the  top  rims  of  the  cylindrical 
cavity.  For  Model  III,  Nodal  Point  193  is  in  the  half¬ 
space,  172  in  the  layer  interface,  and  152,131,111,90,  and 
70  in  the  weathering  layer  as  shown  in  Figure  7c.  The 
wave  arrivals  P^,  S^,  P^,  S^, 


3b  3t  3b  3t  3b  3t 

v2  f2  2'  2  2  2  2  ' 


P2bP2tP2bp2ts2'  PJP2bp2tpjP2t  and  the  -^ultiPles  of 


3b  3t  3b  3t 
2  2  2  2  ' 


P2bP2tP2t’P2ts2  and  are  a11 


identified  on  every  seismogram  of  Figure 


13  a.  and  a0. 


Model  I,  and  II.  In  Figure 
Nodal  Point  193  include  P 


3t 


13  a,,  the  arrivals  at 
_3b  c3f  0  _3t  „3b 

*2  '  b2  '  FJ  2  '  b2  ' 


PfP2tp2bp2tV  PfP3tp2-  PfP3ts2'  SfS3tp2 '  SfSfS2  and 

the  multiples  of  P^P^^'  P3bp3tp2'  P2bp2ts2r  P3bp3ts2 
and  PjP^P^PjP^. 


Nodal  Point  193  is  located  at  a  distance  of  2.9  feet 

from  the  top  rim  of  the  cavity  and  at  a  distance  of  8.6 

feet  from  the  bottom  rim  of  the  cavity.  Since  Nodal  Point 

193  is  located  above  the  cavity,  the  applied  stress  on  the 

top  end-face  of  the  cavity  contributes  to  the  radial 

displacement  at  Nodal  Point  193  in  a  radially  outward 

direction  for  the  arrival  of  P^,  and  the  applied  stress 

on  the  bottom  end-face  of  the  cavity  contributes  to  the 

radial  displacement  at  Nodal  Point  193  in  a  radially 

3b 

inward  direction  for  the  arrival  of  P2  .  In  Figure  13 

a-,  the  high  frequency  components  with  a  period  of -1  milli- 
*  • 

second  are  superimposed  on  the  low  frequency  component  as 
shown  in  the  seismogram  of  Nodal  Point  193.  As  Nodal  Point 


33. 


193  is  located  very  close  to  the  top  of  the  cylindrical 
cavity,  the  location  is  very  sensitive  to  the  motion  on 
the  immediate  top  end  of  the  cavity,  therefore,  the  wave 
PjP^  is  the  type  of  the  wave  traveling  periodically 
from  the  top  rim  of  the  cavity  to  the  axis  of  the  cavity 
with  a  period  of  1  milliseconds.  On  the  trace  of  Nodal 
Point  172  in  Figure  13  a^,  the  high  frequency  arrivals 
with  a  period  around  1.2  milliseconds  are  superimposing 
on  top  of  the  low  frequency  response.  The  wave  P^P^*" 
is  the  type  of  the  wave  traveling  periodically  from  top 
to  bottom  of  the  cavity  along  the  vertical  wall  of  the 
cavity.  The  long  period  wave  is  shown  on  all  the  traces 
in  Figure  13  a,.  This  long  period  waves  is  identified 

2k  2t  ^ 

as  P^  P^  ,  traveling  with  the  P  wave  velocity  in  the  shot 
hole  between  the  two  rims  at  2b  and  2t.  The  Nodal  Point 
152,131,111,90,  and  70  are  located  in  the  first  layer  ( 
weathering  layer) ,  all  the  high  frequency  arrivals  are 
filtered  out  in  the  seismogram  of  Figure  13  a3  and  a^'. 
The  weathering  layer  behaves  as  a  low  pass  filter. 


In  Figure 


13  a^'  and 


i  •  i  _3t  3b 

i  seven  arrivals,  P2  ,  P2 


s3t  s3b 

f  1 


3b  3t  3b  3t 
2  P2  2  2  ' 


and  the  multiples  of  pjP2tp2^PJP2t>P2tpjP2t  are  , 

likewise  the  arrival  of  P^*^  a^-so  is  identified  on  each 
seismogram  in  Figure  13  a- ' .  In  Figure  13  a^',  the 
arrivals  of  p|t,  sf ,  pfp^.  S^S^S,.  pfpfpfpfpfp, 


p2bp2tP  p3t  3b  2t  3b 2t 
P3  3  2’  2  P2  P2  P2  P2  2 


2b  2 1 

and  the  multiples  of  p3  p3  are 


identified  on  the  seismogram  of  Nodal  Point  193;  the  arrivals 

,3t„2t. 


of  P 


3t  s3t  D3tD2tD, 
2  '  b 
2b„2t 


P2 Wp2  ^Pi '  P2  Wp2 "^i '  p3bp3tpi  and  the  multiples 


of  P^P^  are  identified  on  the  trace  of  Nodal  Point  172. 


For  the  rest  of  the  seismogram  the  first  arrival  of  P2  P^ 
is  clearly  identifiable.  In  Figure  13  a^'  and  a^' , 
the  wave  with  a  period  of  about  1.6  milliseconds  as  marked 
by  "A"  on  the  seismogram  of  Nodal  Point  172  in  a2'  is 
identified  as  P  P^P^P  p^fc  on  the  basis  of  the  travel  path 
and  the  travel  time.  The  wave  pjP2  P2  PJP2  t^ie 

of  waves  traveling  periodically  from  bottom  end-face  of 
the  cavity,  through  vertical  face  of  the  cavity,  then  to 
the  top  face  of  the  cavity  with  a  travel  time  of  1.6 
milliseconds.  The  direct  shear  wave  S2  is  clearly 
identified  on  the  seismogram  of  Nodal  Point  193  for  Model 
III,  Figure  13,  a.',  a_',  and  a,'.  The  long  period 
waves  in  Figure  13  a^  and  a^ '  are  P^  P^  p2 »  which 
travels  with  a  P  wave  velocity  of  the  filled  drilling  hole 
periodically  along  the  vertical  face  of  the  hole  between 
the  top  and  bottom  rims  at  2b  and  2t.  The  shape  of  the 
input  forcing  function  is  well  preserved  on  each  seismo¬ 
gram,  in  Figure  13  a^,  a2 ,  a3 ,  a^ ' ,  a2 ' ,  and  a3'  as 
indicated  by  dash  lines. 

.Line  E  > 

Figure  14  shows  the  synthetic  seismograms  of  the 
radial  displacement  Ur  and  the  vertical  displacement  U2 
of  Line  E  at  Nodal  Point  314,335,357,378,400,421,  and  443 
as  shown  in  Figure  7.  Line  E  is  situated  diagonally 
with  respect  to  the  bottom  rim  of  the  cylindrical  cavity. 
Therefore,  the  radial  and  vertical  displacements  shown 
in  Figure  14  a^  and  a^ '  are  identical  to  those 
displacements  in  Figure  13  a^  and  a^ ' ,  except  the 
vertical  displacements  in  Figure  14  a^'  are  in  the 
opposite  direction  and  are  negative  as  shown  in  Figure 
13  a^'.  The  seismograms  as  shown  in  Figure  -14  a2 


and  a, '  are  identical  to  those  as  shown  in  Figure  14 

^  3t 

a-^  and  a^ '  until  the  arrival  of  P2  Pj  in  Model  II,  the 

half-space.  Since  the  nodal  points  of  Line  E  are  located 

in  the  half  space  and  far  away  from  the  layer  interface 

and  the  filled  drilling  hole,  the  effects  of  weathering 

layer  and  the  filled  drilling  hole  on  the  early  response  of 

the  radial  and  vertical  displacement  at  these  nodal  points 

of  Line  E  are  insignificant.  However,  arrivals,  p2tp2» 

p3b_  p3t  P^s  c3tP  p2b_2tp  S3^S  p^bp2t_ 

2  y2'  P2  S2 '  P2  S2'  b2  P2’  P3  P3  2 '  2  fa2 '  3  P3  2' 

P^P^P31^,  P^bP2tP3bS2  and  S3bS3tP2  are  identifiable  in 

Figure  11-14  a3  and  a3 '  . 


Line  F 


Figure  15  shows  the  synthetic  seismograms  of  the 
radial  displacement  Ur  and  the  vertical  displacement  Uz  of 
Line  F  at  Nodal  Point  192,150,108,  and  66. (Figure  7)  For 
Model  III,  Nodal  Point  192  is  in  the  half-space,  and  150 
108,  and  66  in  the  weathering  layer.  In  Figure  15  a. 
and  a2  »  the  arrivals  of  the  waves  p2fc,  P2b,  S2b,  S2b,  P2t>P2t 

and  pjP2^P2tpjP2t  and  bbe  muitiples  P^P^r  and 

PTP^bP~tP7P^t  are  identified.  In  Figure  15  a~,  the 

J  7  ,  „3t  „3t  „2b„ 2t  ^3b^3t  T,3b„3t„2t 

arrivals  of  the  waves  p2  /  S2  '  P2  P2'  P2  P2 '  P2  P2  P2  ' 


P3bP3fc  and  the  multiples  of  P^P^,  P2bP3t,  and  P^P31  are 

also  identified  in  the  seismogram  of  Nodal  Point  192,  which 

is  in  the  half-space.  The  arrival  of  P^,  P^P^,  p?^p?t' 
o^t-^t.-,  _3t_2tc,  n^t^t-.'  _3t_2tc,  _3b_3  tD2  tD  1  . 


P2  P2  Pl'  P2  P2  Sl'  P3  P3  Pl'  P3  P3  Sl'  P2  P2  P2  Pl- and 
the  multiples  of  P^P^,  P^P^,  p2tp2tpi'  p?tp?tsi^  ^nd 

3b  3t  2t  1  2  2  22 

P~  Pt  P,  are  identified  on  the  seismogram  of  Nodal  Point 
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150.  Since  Nodal  Point  150,108,  and  66  are  all  located 
in  the  weathering  layer, which  behaves  as  a  low-pass  filter, 
all  the  high  frequency  responses  are  filtered  out  in  the 
the  seismograms  of  Figure  15  a3  and  a^.  Several  beats 
are  observed  in  the  seismograms  of  the  the  radial  displace¬ 
ment  Ur  in  Figure  15  a^,  a 2,  and  a^. 

The  beat  phenomenon  is  commonly  observed  at  two 
wavetrains  of  slightly  different  frequencies  beat  together 
to  give  the  beats  and  the  amplitude  modulation.  In  the 
present  case,  from  seismogram  at  Nodal  Point  192  in  Figure 
15  a^,  the  period  of  the  envelope  of  the  beats  as  marked 
by  "beats"  in  the  time  interval  of  6.8  milliseconds  and 
11-8  milliseconds  is  measured  to  be  5  millisecon<3s  <  and 
the  period  of  the  amplitude  modulation  on  the  seismogram 
of  Nodal  Point  192  in  the  time  interval  of  8.6  milliseconds 
to  9.3  milliseconds  is  measured  to  be  0.7  milliseconds 
Therefore,  the  frequency  of  the  envelope  of  the  beats  and 
that  of  the  amplitude  modulation  are  found  to  be  about 
196  and  1333  Hz,  respectively.  By  knowing  these  two 
frequencies,  the  periods  of  these  two  wavetrains  are  1.3 
milliseconds  and  1.7  milliseconds.  Therefore,  these  two 
corresponding  waves  are  identified  as  P2  P2  '  travels 

periodically  from  the  bottom  to  the  top  of  the  cavity  along 
the  vertical  face  of  the  cylindrical  cavity  with  a  period 
1.2  milliseconds  and  as  p  P^p^tp  p3t  ,  which  travels 
back  and  forth  periodically  from  the  center  of  the  bottom 
(or  top)  end-face  of  the  cavity,  to  the  rim  ,  then  travels 
along  the  vertical  face  of  the  cavity,  to  the  center  of  the 
top  (or  bottom)  end-face  of  the  cavity  through  the  rim  with 
a  period  of  1.6  milliseconds  These  beat  phenomena  are 
also  clearly  observed  on  the  radial  displacements  of  Nodal 
Point  150,108,  and  66.  A  typical  beating  phenomenon  is 
shown  in  the  seismograms. 


(- 


. 
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In  Figure 


15  a,  '  and  a. 


,3t 


,3t 


P3bP3tS 
P2  P2  b2' 


3t  3b  3t 
*2  P2  2  ' 


, 1 ,  the  arrivals  of  the  waves 
P2tp2bp2ts2-  PJP2bp2tpJP2t' 


and  the  multiples  of 


P2bp2ts2'  P2tp2bp2ts2  and 


p  p-^^p^tp  p3t 

J  2  *2  PJ  2 
nodal  point. 


are  identified  on  each  seismogram  of  the 

3t 

The  arrivals  of  P2  P2  are  also  identified 


on  each  seismogram  as  shown  in  Figure 
15  a-j*/  the  arrivals  of  P2b,  S^, 


15  a 


In  Figure 


P2tp2'  P23ts2'  P2bp2tp2' 


p2bp2ts2'  pfp3tp2'  p3bp3ts2'  and  the 

P^bp3t,  P^P^  are  identified  on  the 

^  ^  ^  ^  3t  2t 

seismogram  of  Nodal  Point  192.  The  arrivals,  P2  P2  P^  , 


multiples  of  P^P^' 


P2  p2  bl' 


S3bs^S  ’ 
b2  b2  bl' 


P3bP3tP|,  P3bP3tS|  and  S^S^PJ  are 


identified  on  the  seismogram  of  Nodal  Point  150.  Nodal 

Point  108  and  66  are  located  in  the  weathering  layer,  only 
3t  2t 

the  arrival  of  P2  P2  P|  is  shown.  In  Figure  15  a^’ 
and  a2 ’ ,  the  wave  in  each  seismogram  with  a  period  very 
close  to  1.6  milliseconds  as  marked  by  "A"  for  Nodal  Point 
192  in  the  time  interval  of  6.8  milliseconds  to  8.4  milli¬ 
seconds,  is  identified  as  PjP^P^PjP^ ,  a  wave  traveling 
periodically  from  the  bottom  end-face  of  the  cavity,  through 
the  vertical  face  of  the  cavity  to  the  top  end-face  of  the 
cavity.  In  Figure  15  a^1,  on  the  seismogram  of  Nodal 
Point  192,  150,  108,  and  66,  the  wave  with  a  period  close 
to  5  milliseconds  as  marked  "B"  on  192  in  the  time  interval 

of  6  milliseconds  and  11  milliseconds  is  identified  as 
2b„2t 


p — p 

P3  P3 


with  a  velocity  of  the  filled  drilling  hole.  In 
Figure  15  a^,  a2 ' »  and  a^',  the  form  of  an  input  step 
pulse  is  well  preserved  in  each  seismogram,  the  amplitudes 
of  the  vertical  displacement  are  larger  than  these  of  the 


radial  displacements.  The  arrivals  of  S 
for  Nodal  Point  192  in  Figure  15  a 
well  identifiable. 


I 

‘1  ' 


3t 
2 
a 


on  the  seismogram 


I 

2  ' 


and  a 


.3 


are 


Line  G 

Figure  16  shows  synthetic  seismograms  of  the  radial 
displacement  Ur  and  the  vertical  displacement  Uz  of  Line  G 
at  Nodal  Point  191,  170,  149,  128,  107,  86,  and  65.  For 
Model  III,  Nodal  Point  191  is  in  the  half-space,  170  in 
the  layer  interface,  and  149,  128,  107,  86,  and  65  in  the 
weathering  layer.  The  locations  of  these  nodal  points  are 
along  the  Line  G  as  shown  in  Figure  7.  The  locations 
of  every  nodal  point  along  Line  G  are  very  close  to  the 
locations  of  nodal  points  along  Line  F.  Line  F  is  in  the 

same  direction  but  with  a  distance  of  0.5  feet  away  from  Line 
G.  The  wave  forms  of  the  radial  displacements  at  Nodal 
Point  191,  149,  107,  and  65  in  Figure  16  a^,  &2>  anc^  a3 
are  very  similar  to  these  of  the  radial  displacement  at 
Nodal  Point  192,  150,  108,  and  66  as  shown  in  Figure  15 
a^  *  a2 1  an^  a3»  except  the  amplitudes  of  the  radial  displace¬ 
ment,  for  Nodal  Point  192,  150,  108,  and  66  are  smaller  than 
these  for  Nodal  Point  191,  149,  107,  and  65.  There  are 
high  frequency  component  responses  observed  at  Nodal  Point 
170,  128,  and  86.  In  examining  the  configuration  of  the 
finite  element  mesh  around  Nodal  Point  149,  170,  and  191, 
as  shown  in  Figure  5,  these  high  frequency  components 
may  be  resulted  from  the  different  arrangements  of  the 
finite  element  mesh  at  Nodal  Point  170  and  Nodal  Point  191 
as  shown  in  the  following  diagram: 
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So  that  Nodal  Point  170  responses  to  a  high  frequence 
radial  displacement  more  efficiently  than  Nodal  Point 
191  or  149.  At  Nodal  Point  170,  it  can  be  viewed  as  if 
two  horizontal  springs  connected  with  two  adjacent  nodal 
points,  contrarily  to  that  at  Nodal  Point  191,  there  are  as 
if  four  diagonal  springs  intersected  at  the  Nodal  Point. 


In  Figure 


16  a. 


and  a. 


the  arrivals  of  the 


waves  P3t,  S3t,  S3b,  P3tP3bP3t,  P2tp2bp2ts2'  and 

PjP^P^PjP^  and  the  multiples  of  P3tP3bP3t,  P^P^P^^1 


3t  3b  3t 
P2  2  P2  b2' 


and  PjP^P^PjP^  are  identified  on  the  seismograms  for 
Nodal  Point  191,  170,  149,  128,  107,  86,  and  65.  In 
Figure  16  a^ ’ ,  for  Model  III,  the  locations  of  Nodal 
Point  191,  170,  149,  128,  107,  86,  and  65  are  in  different 
regions, therefore,  the  arrivals  of  the  waves  are  identified 
each  Nodal  Point  by  each  Nodal  Point.  The  arrivals  of 


waves  P2b, 

g3bg3bp 
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P2bP2tS 
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p3bp2bs 
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p3tp2t 
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P3tP2tS 
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S3tS2tP 
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sfsfs2. 


are  identified  for  Nodal  Point  191.  The  arrivals 


of  p2  ,  p: 

s3bs3tp 

2  2  P3 

P2bP2tP* 
P2  2  1' 


P3bP3tP 
r  2  2  2' 


P3bP3tS 
*2  2  2' 


are  identified  for  Nodal  Point  170, 


_3b  2t 
p2  P2  P3 ' 

The  arrival. 


S2bS2tP' 
b2  b2  *1' 


.3t„2t 


Pi  ,  P 


3t_2t. 


P2bp2tp3' 

p3tp2t  ! 

P3  P3  bl' 


P^hp^S^,  and  S2bS2tS^, 


s3bs2bs  * 

b3  b3  S1 


s3ts2ts 

i>2  b3' 


P3tP2tS 
P3  P3  b3 ' 


S3tS2tP 1 
b3  b3  Pl' 


are  identified  on  the  seismogram  of  Nodal  Point 


149.  The  arrivals  P^P^P.' ,  P3tP2tS.’ ,  P3bP2tP.| ,  P3tP2tP  , 

3t  2t  zzi  i  /.  i  jjj. 

P-.  Si  are  identified  on  the  seismogram  of  128...  Only 
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3  J  2t  , 

the  arrival  of  P2  P2  pi  clearlY  shown  in  the  seismograms 

of  Nodal  Point  107,  86,  and  65. 
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Because  Nodal  Point  191  is  located  right  on  the 
vertical  interface  of  the  filled  drilling  hole,  the 
location  is  very  sensitive  to  the  motion  on  the  top  end 
of  the  cavity.  The  seismogram  of  Nodal  Point  191  in 
Figure  16  a^'  shows  that  the  high  frequency  components 
are  superimposed  on  the  long  period  components.  These 
high  frequency  components  with  a  period  of  about  0.5 
milliseconds  as  marked  by  "A"  on  the  seismogram  of  Nodal 
Point  191  in  a3 '  appear  in  the  time  interval  between  10.9 
milliseconds  and  11.4  milliseconds.  As  the  traveling 
path  of  Pj  is  0.5  foot,  and  the  top  end- face  of  the 
cylindrical  cavity  is  in  contact  with  the  filling  of  the 
drilling  hole  (a  P  velocity  of  1000  ft/sec) ,  the  high 
frequency  components  likely  are  excited  by  P_,  traveling 
periodically  on  the  top  (or  bottom)  end- face  of  the 
cylindrical  cavity.  On  the  other  hand,  the  long  period 
component  in  all  the  seismograms  in  a,'  can  be  identified 

2k  2t  ^ 

as  P3  P3  ,  traveling  periodically  between  the  top  and  bottom 
rims  as  designated  by  2t  and  2b  along  the  vertical  interface 
of  the  drilling  hole  with  a  P  wave  velocity  in  the  drilling 
hole  of  1000  feet/sec  and  a  period  of  5  milliseconds. 

Nodal  Point  149,  128,  107,  86,  and  65  are  located  in  the 
weathering  layer,  therefore,  most  of  the  high  frequency 
response  are  apparently  filtered  out. 


Line  H 


Figure  17  shows  the  responses  of  both  the  radial 
component  and  the  vertical  component  at  Nodal  Point  190, 
140,  106,  and  64  along  Line  H.  All  the  nodal  points 
along  this  line  are  on  the  top  region  of  the  cylindrical 
cavity  source.  For  Model  III,  these  points  are  located 
in  the  filled  drilling  hole.  Similar  to  that  of  Line  G 
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and  Line  F,  a  beating  phenomenon  is  again  observed 
in  the  time  history  of  the  radial  displacement  at 
these  four  nodal  points  along  Line  H  for  Model  I  and 
II.  In  the  present  case,  these  beats  have  much 
smaller  amplitude  modulations  .than  these  of 
Lines  F,  and  G.  As  nodal  points  along  Line  H  is  located 
close  to  the  axis  of  the  cylindrical  cavity,  therefore, 
the  radial  motion  of  the  top  end-face  of  the  cavity 
contributes  to  the  vertical  displacement  significantly. 

The  high  frequency  response  due  to  the  radial  motion  of  the 
end-faces  of  the  cavity  is  clearly  observed  on  the  synthetic 
seismograms.  The  wave,  P  P^P^P  P^fc,  is  clearly 


identified  in  the  seismogram  of  the  vertical  displacement 
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U2,  which  are  very  similar  to  that  of  PjP^  P^  PjP2  *-n 
Figure  15,  and  Figure  16.  For  Model  III,  Figure 
17  a-j'r  because  the  nodal  points  are  located  in  the 
drilling  hole,  which  is  filled  with  unconsolidated  material, 
the  amplitudes  of  both  the  radial  and  the  vertical 
displacements  are  very  much  amplified,  for  instance  an 
extremely  large  amplitude  wave  is  observed  at  Nodal  Point  190. 
The  wave  is  attenuated  rapidly  as  the  observation  point 
moves  away  from  the  cavity  (see  both  in  Figure  17  a. 


and  a, ' ) . 

J_2b_2t 
waves  P3  P3 


In  Figure 


17  a,,  at  Nodal  Point  148,  the 


2  v  2  ^  J 

and  P„  pf  apparently  travel  back  and  forth 


periodically  along  the  vertical  interface  between  the  two 

rims  as  designated  by  2t  and  2b.  At  Nodal  Point  148,  106 

and  64,  the  most  pronounced  arrivals  are  indicated  on  the 

3t  2t  2t 

seismograms.  The  arrival  of  P3  P3  P^  with  a  period  of 
about  5  milliseconds  is  indicated  on  each  seismogram  of 
Nodal  Point  148,  106,  and  64. 
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Line  I ,  J ,  K 


Figures  18,  19,  and  20  give  synthetic  seismograms 
of  the  radial  and  the  vertical  displacements  at  the  nodal 
points  along  Line  I,  J,  and  K.  These  three  lines  are 
located  in  the  region  bellow  the  bottom  of  the  cavity. 
Geometrically  the  locations  of  these  nodal  points  along 
these  three  lines  are  symmetric  to  those  of  nodal  points 
along  Line  F,  G,  and  H.  For  the  whole-space  case,  the 
symmetrical  responses  of  the  radial  and  vertical  displace¬ 
ments  are  displayed  in  a^  and  a^  of  Figure  18,  19,  and 
20,  and  can  be  compared  with  the  results  of  a.^  and  a^  of 
Figure  15,  16,  and  17  to  be  exactly  identical.  For 
the  half-space  case  as  compared  with  the  whole-space  case, 

the  only  difference  is  the  arrival  of  additional  reflected 
3t 

wave  P 2  P2  from  the  surface  of  the  half-space  as  indicated 
on  each  seismogram  for  each  nodal  point  as  shown  in  a2  and 
a2  of  Figure  18,  19,  and  20.  For  Model  III,  the  layered 
half-space  case,  the  synthetic  seismograms  in  general  are 
quite  similar  in  appearance  as  these  for  the  whole-space  case, 
except  the  additional  diffracted  arrivals. 

r 

All  the  above  interpretations  are  limited  to  the 
arrivals  in  the  time  interval  of  0  to  15  milliseconds. 
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7.  CONCLUSION 


The  present  study  has  demonstrated  the  possibility  of 
obtaining  the  numerical  solutions  for  a  finite  length  energy 
source  within  a  half-space,  and  a  layered  half-space  by  means 
of  the  finite  element  method.  The  significant  results 
regarding  the  displacement  field  due  to  a  finite  cylindrical 
cavity  source  in  a  whole-space,  a  half-space,  and  a  layered 
half-space,  may  be  summerized  as  follows: 

1)  .  For  the  observation  points  inside  the  medium,  the 
reflected  waves  from  the  surface  boundary  of  the  half-space 
for  the  velocity  of  the  medium  assumed  does  not  seem  to  have 
significant  effects  on  the  amplitude  of  the  wave. 

2)  .  The  amplitudes  of  both  the  vertical  and  radial 
displacements  at  the  observation  point  on  the  surface  of  the 
half-space  model  are  much  larger  than  these  at  these  corres¬ 
ponding  observation  points  in  the  whole-space. 

3)  .  The  wave  forms  of  the  transient  response  strongly 
depend  on  the  dimensions  of  the  finite  cylindrical  cavity,  i.e. 
the  radius  of  the  cavity,  and  the  length  of  the  cavity. 

4)  .  The  type  of  P  wave  traveling  periodically  from  the 
top  rim  to  the  bottom  rim  of  the  cavity  along  the  vertical 
interface  of  the  cylindrical  cavity  and  the  type  of  P  wave 
traveling  periodically  from  the  bottom  (or  the  top)  face  of 
the  cavity  with  a  motion  described  by  the  Bessel  function 
order  zero,  through  the  vertical  face  of  the  cavity,  to  the 
top  (or  the  bottom)  face  of  the  cylindrical  cavity  also-  with 
a  motion  described  by  the  Bessel  function  order  zero  are 
dominanted  in  all  three  models. 


5)  .  There  are  no  vertical  displacements  L'z  at  any  point 
along  the  middle  plane,  which  is  symmetrical  with  respect  to 
the  two  end  faces  of  the  cylindrical  cavity,  in  Model  I  -  a 
whole-space . 

6) .  For  Model  III  -  a  layer  half-space,  the  long  period 
wave  is  clearly  identifiable,  traveling  periodically  along 
the  vertical  face  of  the  drilling  hole  between  the  two  rims, 
the  rim  of  the  contact  at  the  layer-interface  and  the  top 

rim  of  the  cylindrical  cavity. 

7)  .  All  the  high  frequency  arrivals  are  filtered  out 
in  the  seismograms  at  those  observation  points  located  in  the 
weathering  layer.  The  weathering  layer  behaves  as  a  low 
pass  filter. 

8 ) .  The  shape  of  the  input  forcing  function  is  will 
preserved  on  each  seismogram  of  the  vertical  displacement 
Uz  but  not  the  radial  idsplacement  Ur. 

i 

.  9) .  The  beat  phenomenon  is  clearly  observed  on  the  radial 

^displacements  of  these  observation  points  located  in  the  upper 
region  of  the  cavity  source. 

10)  .  The  amplitudes  of  both  the  radial  and  vertical  dis¬ 
placements  are  extremely  large  at  these  observation  points  close 
to  the  cavity  in  the  filled  drilling  hole,  but  are  attenuated 
rapidly  as  the  observation  point  moves  away  from  the  cavity. 

11)  .  Shear  wave  can  be  generated  from  the  finite  length 
cavity  source  with  only  normal  stresses  applied  on  the" 
surfaces  of  the  cavity.  The  arrival  of  shear  waves  are 
clearly  identifiable  at  these  observation  points  in  the  region 
along  the  diagonal  line  of  the  cylindrical  cavity. 
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ABSTRACT 


Three  dimensional  model  experiments  are  performed  to  study  the 
acoustic-elastic  wave  scattering  and  diffraction  by  a  vertical  elastic 
cylinder  immersed  in  a  fluid  due  to  a  trainsient  acoustic  point  source. 

The  direct,  transmitted  and  diffracted  waves  in  the  volume  exterior  to 
the  cylinder  are  clearly  observed.  Particularly,  the  diffracted  longi¬ 
tudinal,  transvers  •,  Rayleigh  and  acoustic  (creeping)  waves  which  propa- 
age  near  the  cylindrical  interface  are  observed  together  with  the 
waves  transmitted  through  the  cylinder.  The  moveout  of  the  diffracted 
waves  for  non- symmetric  source-receiver  configurations  is  shown. 

Splitting  of  each  diffraction  at  non- symmetric  receiver  positions  is 
verified  and  explained  as  a  result  of  diffracted  energy  propagating  in 
both  directions  around  the  cylinder.  Helical  diffracted  waves  propa¬ 
gating  up  and  down  the  cylinder,  as  predicted  by  theory,  are  verified 
for  the  first  tine. 

Three  dimensional  ray  path  and  wavefront  equations  are  given  i 01 
all  of  the  waves  encountered  in  this  experiment.  The  amplitude  spectra 
maxima  are  shown  to  be  related  to  the  relative  arrival  times  between 
the  various  events  observed  on  the  data  trace.  Envelope  dectection  and 


homomorphic  filtering  ire  applied  and  shown  to  provide  information 


or.  the  arrival  times  and  cn  the  temporal  duration  of  the  observed 
transmitted  and  dj  lira,  red  event:;.. 
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INTRODUCTION 


Bodies  possessing  a  cylindrical  shape  are  of  considerable  impor¬ 
tance  in  the  theory  of  scattering  and  diffraction  of  acoustic  and 
elastic  waves.  Many  man-made  structures,  such  as,  mine  shafts  and 
drifts,  fiber  reinforced  materials,  underground  and  ocean-bottom 
pipes,  cast  metals  containing  inclusions  or  imperfections  and  missile 
silos,  as  well  as  such  natural  formations  like  salt  domes,  mineral 
deposits,  igneous  intrusives,  caves,  sink  holes,  seamounts  and  coal 
gasification  cavities  may  be  approximated  geometrically  by  a  cylindri¬ 
cal  shape.  Since  the  analytical  solution  to  the  three  dimensional 
scattering  and  diffraction  of  an  elastic  cylinder  embedded  in  either 
an  elastic  or  an  acoustic  medium  due  to  a  transient  point  source  of 
energy  is  intractable,  the  alternative  method  of  an  experimental  model 
study  provides  a  means  of  obtaining  a  physical  understanding  of  the 
problem.  Moreover,  an  experimental  study  may  provide  a  means  to  ver¬ 
ify  that  the  waves  predicted  by  two  dimensional  theory  do  exist  in 
both  two  and  three  dimensions. 

The  basic  physics  of  scattering  and  diffration  at  various  discon¬ 
tinuities  is  dicsussed  by  Keller(1958,  1962),  Keller  and  Karal(1964) 
and  Kouyoumjian  (1975).  Excellent  reviews  of  the  scattering  and  dif¬ 
fraction  of  fluid-loaded  structures  are  provided  by  Uberall  and 
Hwang(1976),  Uberall(l973)  and  Neubauer(1973) .  These  reviews  are  lim¬ 
ited  to  discussing  the  analytical  solutions  to  simple  acoustic  and 
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elastic  problems  by  the  method  of  integral  transforms  and  they  discuss 
the  use  of  surface  waves  propagating  on  plane  and  curved  boundaries. 
The  principal  analytic  solutions  which  have  been  obtained,  only  begin¬ 
ning  to  scratch  the  surface  of  acoustic  and  elastic  scattering  and 
diffraction  by  cylindrical  and  spherical  bodies,  are  for  sound  pulses 
diffracted  by  a  circular  cylindrical  cavity(Chen  and  Pao,  1977),  for 
sound  pulses  diffracted  by  a  rigid  and  a  soft  cylinder  embedded  in  a 
solid (Gilbert , 1960  and  Gilbert  and  Knopoff,  1959),  for  three  dimen¬ 
sional  acoustic  wave  scattering  and  diffraction  by  an  open-ended  ver¬ 
tical  soft  cylinder  in  a  half-space(Teng,  Kuo  and  Gong, 1975)  and  for 
the  scattering  and  diffraction  of  an  infinite  cylindrical  cavity  in  an 
elastic  medium  due  to  an  impulsive  elastic  P  wave  source(Hwang  and 
Kuo,  unpublished).  These  solutions  provide  a  basis  for  the  under¬ 
standing  of  the  manner  in  which  incident  impulsive  energy  is  scattered 
and  diffracted  by  cylinders,  but  they  fall  short  of  the  total  problem 

of  showing  the  interaction  of  the  transmitted  waves  through  an  elastic 

\ 

cylinder  with  the  diffracted  waves  around  the  cylinder.  Failure 
occurs  because  the  boundary  conditions  used  prohibit  energy  from  being 
refracted  into  the  cylinder.  Incorporation  of  transmission  into  the 
problem  leaves  a  characteristic  equation  to  solve  whose  roots  are 
extremely  difficult  to  extract.  Except  for  the  solutions  by  Teng  and 
Kuo(1975)  and  Hwang  and  Kuo (unpublished) ,  the  previously  mentioned 
solutions  are  for  two  dimensional  problems. 
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Experimental  model  studies  have  provided  the  verification  of  the 
existence  of  circumferential  diffracted  creeping  waves  and  of  dif¬ 
fracted  Rayleigh  waves  (Barnard  and  McKinney,  1961,  Faran,  1951,  Bun- 
ney,  Goodman  and  Marshall,  1969,  Harbold  and  Steinberg,  1969,  and  Neu- 
bauer,  1968).  These  experiments  were  specifically  designed  to  study 
individual  diffractions  by  using  a  bistatic  source  and  receiver 
configuration  with  the  transducers  oriented  to  the  cylindrical  surface 
at  the  critical  angle  at  which  the  diffracted  wave  to  be  verified  is 
excited.  Through  this  orientation  the  desired  diffraction  event  is 
enhanced  while  all  of  the  other  arriving  events  are  less  enhanced  in 
order  that  they  will  not  interfere  and  confuse  the  verification  of  the 
diffraction.  Except  for  the  work  of  Steinberg(1969) ,  all  of  the 
experiments  were  performed  for  two  dimensional  cylinders, i.e.  the 
source  and  receiver  were  kept  coplanar.  Another  drawback  of  the  pre¬ 
vious  experiments  designed  to  study  diffractions  on  cylinders  is  that 
long  monochromatic  tonebursts  were  used  as  sources  which  were  long 
compared  to  the  characteristic  dimension  of  the  cylindrical  scatterer. 
Sachse(1974) ,  Sachse  and  Chian(1974)  and  Pao  and  Sachse(1974)  used  a 
short  transient  pulse  to  study  backscattered  and  diffracted  ultrasonic 
pulses  from  elastic  and  fluid  cylindrical  inclusions  in  a  solid  med¬ 
ium.  Verification  of  a  circumferential  diffracted  S  wave  propagating 
at  a  fluid-solid  interface  for  the  case  of  a  fluid  inclusion  in  a 
solid  medium  was  obtained. 
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This  thesis  addresses  the  acoustic-elastic  scattering  and  diffrac¬ 
tion  of  an  infinite  elastic  cylinder  immersed  in  a  fluid  medium  due  to 
a  simulated  transeint  point  source. 


ACOUSTIC-ELASTIC  EXPERIMENTAL  MODEL 


Model  System 

The  experiments  are  performed  in  an  acoustic  model  tank  which  mea¬ 
sures  24'  long  by  18'  wide  by  12'  deep  (Figure  1).  By  suspending  the 
source,  the  receiver  and  the  elastic  cylindrical  model  in  the  center  3 
sq.  ft.  of  the  tank,  the  sides  and  the  bottom  are  far  removed  to  pre¬ 
vent  the  recording  of  unwanted  reflections.  In  the  present  experi¬ 
ment,  using  a  rectangular  coordinate  system,  the  source  and  the 
receiver  can  be  moved  in  both  directions  along  a  system  of  rails  using 
the  position  of  the  cylinder  as  the  origin.  For  the  cylindrical  sys¬ 
tem,  the  source  transducer  is  used  as  the  origin  and  the  receiver  is 
rotated  about  the  source  at  a  radius  r.  In  both  systems  the  vertical 
movements  are  achieved  by  raising  or  lowering  the  transducers  attached 
to  the  rods.  Both  the  source  and  the  receiver  can  be  rotated  verti¬ 
cally  and  horizontally  about  their  fixed  positions  in  space  to  simu¬ 
late  a  point  source  or  a  line  source.  Through  the  rotation  of  the 
source  in  a  horizontal  plane  in  degree  increments  and  summing  the 
recorded  data  traces  for  each  increment,  a  line  source  can  be  simu¬ 
lated  because  a  narrow  beam  transducer  is  used  as  the  source.  If,  for 
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each  horizontal  increment,  the  source  is  rotated  vertically  in  equal 
increments  and  these  data  are  summed,  a  point  source  is  simulated. 

The  source  and  the  receiver  are  suspended  at  a  depth  of  one  foot  or 
greater  below  the  free  surface  to  eliminate  unwanted  reflections  from 
the  surface.  A  cylinder  of  4.5  ft.  in  length  is  used.  The  source  and 
the  receiver  are  kept  far  enough  from  each  end  of  the  cylinder  so  that 
diffrations  from  the  ends  do  not  interfere  with  the  arrivals  of  pri¬ 
mary  interest. 

The  source  and  the  receiver  used  in  the  experiments  are  immersion 
transducers  manufactured  by  Panametrics,  Inc..  These  transducers  mea¬ 
sure  1.25  in.  in  length  and  0.625  in.  in  diameter.  The  element  used 
is  a  0.5  in.  diameter  PZT  (Lead-Zirconium-Titanate)  ceramic  which  is 
backed  by  a  specially  designed  tungsten-doped  epoxy  resin.  Figure  2b 
displays  a  typical  source  wavelet  for  this  type  of  transducer  in  The 
pulse  mode  along  with  its  amplitude  and  phase  spectra  as  shown  in  Fig¬ 
ure  2a.  The  pulse  emitted  from  the  transducer  used  as  the  source  is 
short,  approximately  2  psec.  duration,  and  the  radiation  patterns 
are  about  10  db.  down  at  10°  from  the  transducer  axis  creating  a  nar¬ 
row  directed  beam  of  energy  (Figure  3).  The  lack  of  significant  ring¬ 
ing,  as  seen  in  the  wavelet  in  Figure  2b,  when  the  transducer  is 
pulsed  by  a  narrow,  1  psec.,  high  voltage(approximately  80  volt) 
pulse  is  obtained  by  highly  damping  the  transducer  using  a  tungsten- 
doped  epoxy  resin  to  back  the  PZT  ceramic  disk.  This  resin  is 


a.J  Receiver  Radiation  Pattern 


Figure  3 


designed  to  match  the  impedance  of  the  PZT  ceramic  allowing  100%  of 
the  energy  to  be  transmitted  at  the  back  of  the  ceramic.  Since  the 
epoxy  has  a  well  determined  velocity,  by  doping  the  epoxy  with  tung¬ 
sten  dust  the  density  of  the  resin  can  be  increased  enough  that  the 
impedance  of  the  resin  matches  the  impedance  of  the  PZT  ceramic.  How¬ 
ever,  the  epoxy  has  a  low  Q  which  completely  attenuates  the  signal 
transmitted  at  the  back  of  the  ceramic  before  it  reaches  the  ceramic 
following  reflection  from  the  rear  casing.  These  transducers  were 
chosen  specifically  because  they  emitted  a  short  almost  symmetric 
pulse  when  used  in  the  pulse  mode.  This  type  of  pulse  is  desired 
because  it  allows  greater  resolution  between  arrivals. 

Data  Acquisition  System 

The  experimental  data  acquisition  system  is  shown  in  Figure  4.  It 
uses  an  HP-214A  pulse  generator  to  pulse  a  Panametrics  immersion 
transducer  with  a  narrow,  1  psec.,  high  voltage  pulse  to  provide  a 
narrow,  2  psec,  source  wavelet  (Figure  2b).  The  received  signal, 
which  is  scattered  and  diffracted  from  the  cylinder  model,  is  first 
received  by  a  Panametrics  immersion  transducer,  then  it  is  pre-ampli- 
fied,  then  amplified  20  db.  by  an  HP-450A  amplifier.  It  is  then  fil¬ 
tered  by  a  100  KHZ.  high  pass  RC  Rockland  filter  before  it  enters  the 
Tektronix  sampling  unit  for  reduction  to  an  analog  signal  compatible 
with  the  8  bit  A/D  converter.  Since  the  center  frequency  of  the 
transducer  pulse  is  in  the  megahertz  range (see  Figure  2b),  a  Tektronix 
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585  oscilloscope  with  a  1S1  sampling  unit  plug-in  is  used  to  slow  the 
signal  down  to  a  recordable  level.  A  replica  of  the  input  signal  to 
the  sampling  unit  is  output  as  a  lower  frequency  analog  signal  which 
is  digitized. 

The  analog  output  signal  from  the  sampling  unit  is  input  into  an  8 
bit  analog-to-digital  (A/D)  converter  which  has  a  memory  of  1024  8-bit 
words.  Since  the  desired  signal  requires  more  than  1024  samples  to 
record  the  entire  signal,  the  A/D  conveter  has  a  built-in  option  which 
allows  the  recording  of  successive  records  of  1024  samples  by  counting 
(N-l)21024  samples  before  beginning  to  digitize  the  signal  (N  =  the 
number  of  the  data  group  which  is  desired).  After  each  data  group  is 
written  into  the  memory  of  the  A/D  converter,  it  is  read  out  onto  an 
8-bit  paper  tape  which  is  then  read  into  a  PDP  11/50  computer  for  the 
appropriate  pre-processing  prior  to  analysis  on  an  IBM  360-91/75  com¬ 
puter.  The  pre-processing  performed  on  the  PDP  11/50  consists  of 
merging  and  editing  the  appropriate  data  groups  together  to  form  one 
data  record  for  each  receiver  position  recorded.  If  the  data  were 
clipped,  this  record  would  be  read  to  paper  tape  and  the  amplification 
factors  would  be  adjusted  to  record  an  unclipped  record.  Once 
recorded  both  clipped  and  unclipped  records  are  read  into  the  computer 
and  corrected  using  a  scaling  program  to  obtain  the  final  unclipped 
record  of  the  data  trace  for  that  receiver  position. 
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Model 


The  model  used  is  a  2.0  inch,  4.5  feet  long,  stainless  steel  cylin¬ 
der  which  has  a  longitudinal  wave  velocity  of  18996  ft/sec  and  a  tran¬ 
sverse  wave  velocity  of  10170  ft/ sec.  Three  dimensional  scattered  and 
diffracted  fields  of  a  simulated  acoustic  point  source  due  to  a  cylin¬ 
der  are  studied:  (1)  the  2D  coplanar  case;  (2)  the  3D  non-coplanar 
case.  Receivers  are  placed  in  the  same  horizontal  plane  as  the  source 
for  Case  (1)  and  for  Case  (2)  they  are  placed  in  a  horizontal  plane 
different  from  the  source.  The  non-coplanar  case  is  particularly 
designed  to  study  the  helical  waves.  The  coplanar  case  simulates  the 
scattering  and  diffraction  from  an  infinite  elastic  cylinder  immersed 
in  a  fluid  due  to  an  acoustic  line  source  parallel  to  the  cylinder 
axis. 

The  geometrical  configuration  of  the  source  and  the  receivers  for 
both  the  coplanar  and  the  non-coplanar  cases  are  shown  in  Figures  5 
and  6.  The  same  azimuthal  receiver  locations  are  used  throughout  the 
experiment.  Only  the  horizontal  plane  of  the  receivers  is  altered. 
In  order  to  obtain  data  in  the  illuminated  zone  the  source  and  the 
receiver  locations  were  changed  to  6Rq  as  shown  in  Figure  6. 

The  configurations  of  the  source  and  the  receivers  in  this  study 
are  divided  into  four  regions.  Region  I  consists  of  receivers  at  azi¬ 
muths  0°  to  30°  on  a  circular  arc  of  3.5Rq  about  the  cylinder  for  a 
source  located  3.5Rq  from  the  cylinder  axis.  These  receivers  are 


source 


Figure  5 
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Source  and  receiver  configuration  for 
regions  I,  II,  III. 

GSB  =  Geometrical  Shadow  Boundary 
SRB  =  S  Reflection  Boundary 
PRB  =  P  Reflection  Boundary 
SdB  *  Diffracted  S  Wave  Boundary 
PhB  =  Diffracted  P  Wave  Boundary 


Figure  6:  Source  and  receiver  configuration  for 
1  region  IV. 

GSB  =  Geometrical  Shadow  Boundary 
SRB  =  S  Reflection  Boundary 
PRB  *  P  Reflection  Boundary 
S^B  =  Diffracted  S  Wave  Boundary 
PdB  =  Diffracted  P  Wave  Boundary 
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all  located  in  the  geometrical  shadow  zone  beyond  the  geometrical  sha¬ 
dow  boundary(GSB) .  Regions  II,  III  and  IV  are  all  within  the  illumi¬ 
nated  zone.  Region  II  has  the  source  located  at  3.5Rq  with  the 
receivers  also  at  3.5Rq  at  azimuths  35°  to  55°.  This  region  is  bey¬ 
ond  the  S  Reflection  Boundary (SRB)  where  reflected  waves  are  not  dec- 
tected  but  diffracted  waves  are  detected.  Region  III  has  the  source 
located  at  3.5Rq  while  the  receivers  are  on  a  4.5Rq  circular  arc 
at  azimuths  75°  to  105°.  This  region  is  within  the  S  and  P  Reflection 
Boundaries (SRB  and  PRB)  where  critically  reflected  and  diffracted 
waves  propagate.  For  Region  IV,  as  shown  in  Figure  6,  the  source  is 
located  at  6Rq  and  the  receivers  are  placed  at  azimuths  120°  to  160° 
6Rq  from  the  cylinder.  This  region  is  within  the  boundaries  for  P 
and  S  wave  diffractions  where  only  reflected  and  transmitted  waves 
propagate.  Figures  7  display  the  interpreted  data  for  the  Region  I, 
the  shadow  zone(SZ).  Figures  8,9  and  10  display  the  interpreted  data 
for  Regions  II,  III  and  IV  located  in  the  various  parts  of  the  illumi¬ 
nated  zone(IZ). 

ANALYSIS,  INTERPRETATION  AND  DISCUSSION 
Geometry  of  Raypaths (Wavefront  Parametric  Equations) 

Definition  and  notation  of  various  types  of  waves 

The  following  symbols  are  used  to  define  the  various  types  of 
waves  in  this  thesis.  K  is  used  to  denote  the  acoustic  waves  propa¬ 
gating  in  the  acoustic  medium  outside  the  elastic  cylinder.  P  and  S 
are  used  to  denote  the  longitudinal  and  transverse  waves  propagating 
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Bottom  3  traces  pre  2D  data 
Top  3  traces  are  3D  data 


y 

inside  the  elastic  cylinder.  and  denote  the  diffracted 

longitudinal  and  transverse  waves  which  propagate  at  the  cylinder 
interface.  and  R^  denote  the  diffracted  acoustic  wave  and 
the  diffracted  Rayleigh  wave  , respectively. 

Combinations  of  the  symbols  K,  P,  S,  P^,  S^,  R^  and 
are  used  to  describe  the  ray  paths  of  the  waves  interpreted  in  this 
study.  The  direct  acoustic  wave  which  has  no  interaction  with  the 
cylinder  model  is  denoted  by  K.  KK^  denotes  an  incident  acoustic 
wave,  K,  reflected  from  the  cylinder.  KPK  denotes  an  incident 
acoustic  wave,  K,  refracted  into  the  elastic  cylinder  as  a  longitu¬ 
dinal  wave,  P,  and  refracted  into  the  acoustic  medium  as  an  acoustic 
wave,  K.  Analogous  to  KPK  is  KSK  ,  the  difference  being  the  wave 
travels  as  a  transverse  wave,  S,  inside  the  elastic  cylinder.  Waves 
which  are  multiply  reflected  a  P  waves  inside  the  elastic  cylinder 
are  denoted  by  KnPK,  where  n  equal  the  number  of  P  wave  ray  paths. 
Similarly,  KmSK  denote  multiply  internally  reflected  S  waves  with  m 
being  the  number  of  S  wave  ray  paths,  points  inside  the  elastic 
cylinder.  Diffracted  waves  are  denoted  as  follows:  KP^K  is  the 

X 

diffracted  P  wave;  KS^K  is  the  diffracted  S  wave;  KRdK  is  the 
diffracted  Rayleigh  wave;  KK^K  is  the  diffracted  acoustic  wave. 
Full  circumnavigations  of  the  cylinder  by  the  various  diffracted 
waves  are  denoted  by  KP^K  and  similarly  for  the  other  diffrac- 
teins.  For  non-symmetric  source-receiver  configurations  ,  as  shown 
in  Figures  5  and  6  the  diffractions  propagate  in  both  directions. 
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The  front  diffraction(  the  shortest  path)  is  denoted  by  KP^K^ 
and  the  back  diffraction(the  longest  path)  is  denoted  by  KP^K*. 
the  same  symbolism  is  used  for  all  of  the  front  and  back  diffracted 
waves . 

Direct  Waves 

The  direct  waves,  K,are  only  observed  in  the  illuminated  zone, 
i.e.,  at  receivers  in  Regions  II,  III  and  IV.  In  Region  II,  K 
arrives  at  approximately  112  psec  at  35°  coupled  with  KR^jK** 

(see  Figure  8).  At  40°,  K  arrives  coupled  with  K3PK  and  again  at 
55°  coupled  with  KPSK  at  approximately  105.5  psec.  At  45°,  K  is 
not  visible.  There  is  a  data  void  from  107  to  108.25  psec.  Des¬ 
tructive  interference  as  well  as  low  source  strength  account  for  the 
loss  of  data  where  K  should  appear.  Unlike  Region  II,  K  arrives 
among  the  first  arrivals  in  Region  III(see  Figure  8).  It  is  coupled 
with  KPK,  K2PK  and  KP^K^  arriving  at  approximately  109.25  psec 
at  75°  for  the  coplanar  case.  Region  IV  is  devoid  of  K  arrivals  due 
to  the  source  radiation  pattern  and  the  geometrical  configuration 
used. 
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Reflected  Waves 


Region  II  receives  no  reflected  events,  KK^ ,  because  this 
region  is  beyond  the  critical  reflection  boundaries  for  P  and  S 
waves.  Critically  reflected  P  and  S  waves  are  excited  beyond  the  P 
and  S  diffracted  wave  boundaries.  The  most  prominent  KK^  arrivals 
appear  at  receivers  in  Region  IV.  They  are  the  first  arrivals  and 
their  amplitude  is  much  greater  than  the  diffracted  and  refracted 
wave  arrivals.  In  order  to  show  all  of  the  arrivals,  the  KK 

r 

events  are  displayed  clipped. 

The  wavefronts  for  the  reflected  waves,  KK^,  can  be  described 
parametrically.  According  to  Huygens'  Principle  the  surface  of  the 
cylinder  will  have  and  an  infinite  number  of  secondary  sources  when 
the  incident  energy  from  the  source  irradiates  it.  Let 
Q(£,H,C)  be  the  coordinate  point  of  the  secondary  source  on 
the  surface  of  the  cylinder.  Let  P(x,y,z)  represent  a  point  on  the 
advancing  wavefront  of  the  reflected  wave  from  the  secondary  sources 
on  the  cylinder's  surface.  The  length  of  the  incident  ray  from  the 
source,  S(0,0,0),  to  the  cylinder  is  given  parametrically  as 

R^  =  (YQ*  -  2aYgCos0  +  a2sec*$)^  (1) 
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where  Yg  is  the  y  coordinate  of  the  offset  cylindrical  axis  and  a 
is  the  radius  of  the  cylinder.  The  parameters  used  to  describe  the 

wavefront  are  0  and  $(see  Figure  11).  Parametrically,  the  point 
Q(4,n,C)oa  the  surface  of  the  cylinder  is  given  as 

£  =  asin6  (2a) 


H  =  YQ-  acos0  (2b) 

£  =  atan$  .  (2c) 

From  the  surface  of  the  cylinder,  the  general  equation  for  the 
expanding  wavefront  in  three  dimensions  from  Q(£,n>0  is 

(*  -  £)2  +  (y  -  n)2  +  (*  -  C)2  =  h2  (3) 

where  R  is  the  radius  of  the  expanding  wavefront  at  any  point 
P(x,y,z).  In  order  to  describe  the  surface  of  the  wavefront  due  to 
all  of  the  secondary  sources  at  a  given  point  Pg(Xg,yg,Zg) , 
equation  (3)  is  differentiated  with  respect  to  6  and  0  and  then 
these  equations  are  solved  along  with  equation  (3)  for  the  coordi¬ 
nates  x,  y,  and  z.  This  solution  provides  the  coordinates  of  the 
point  P0(x0>y)>2())  on  the  wavefront  of  the  reflected 
waves, KK^.  Differentiating  equation(3)  we  obtain 
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a(x^  -  £)cos0  +  a(yQ  -  q)sin0  =  -R3R/30  (4a) 

a(zQ  -  £)sec2$  =  -R3R/3*  .  (4b) 

The  radius  of  the  expanding  wavefront  for  KKf  is 

Rj  =  c2t  -  R.  (5) 

where  t  is  the  time  variable,  c  is  the  acoustic  velocity  of  the  fluid 
surrounding  the  cylinder  and  R^^  is  the  length  of  the  ray,  or  the 
wavefront  radius,  from  the  source  to  the  cylinder  as  given  in  equation 
(1).  By  substituting  equations  (1),  (2)  and  (5)  into  equations  (3) 
and  (4)  and  solving  for  X,  y  and  z,  we  obtain  the  following  solution 
for  a  point  p0^x0,y0,20^  °n  tbe  wavc^ront 

Xq  =  4  +  *YgSin0cos8  ±  sin0(l  -  Rr  2 

YQ2sin20  -  a2R/"  2tan2<)>)^)  (6a) 


yQ  =  q  +  R1(Ri"1Yosin20  T  cos0(l  -  R-  2 


Yg2sin26  -  a2Rr  2tan2(}»)Ji) 


(6b) 
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zQ  =  C  +  •RjRi’1tan^ 


Equations  (6)  describe  the  three  dimensional  wavefront  of  the  waves 
reflected  from  the  cylinder.  If  0  =  0  equations  (6)  reduce  to  the 
two  dimensional  wavefront  and  ray  path  equations  as  follow: 


Xq  =  4  +  RjCRj  *YgSin8cos6  ±  sin0(l  -  R~  2 


Y02sin28)^ 


yQ  =  n  +  R1(Ri"1Y()sin2e  cos6(l  -  R/“ 2YQ2sin29)^) 


20  *  ° 


R^.  =  (YQ2  +  a2  -  2aYgCOS0) 


Refracted  Waves 


Refracted  waves  are  not  recorded  in  this  study  because  the 
receivers  are  restricted  to  recording  in  the  fluid  exterior  to  the 
cylinder.  The  parametric  wavefront  equations  are  given  for  com¬ 
pleteness.  The  radius  of  the  refracted  waves  expanding  wavefront  is 
given  as  follows 
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Yp  =  v/c. 


(8) 


R2  =  Yp(ct  -  R.) 


where  v  is  the  velocity  of  the  refracted  wave.  For  the  refracted  P 
wave,  KP,  v  =  a  and  for  the  refracted  S  wave,  KS  ,  v  =  p.  The 
length  of  the  incident  ray,  R^,  is  given  by  equation  (1)  and  the 
point  Q(£,n>0  on  the  cylinder's  surface  is  qiven  by  equations 
(2).  Substituting  equations  (1),  (2)  and  (8)  into  equations  (3)  and 
(4)  we  obtain  the  parametric  equations  for  the  wavefront  of  the 
refracted  wave  (P  or  S)  at  Pq  as  follows: 


xQ  =  £  +  R2(Ri~1Y()sin0cose  ±  sin0(l  -  Yp2 


Ri  2^QZs*n2®  “  Yp2*2®^"  2tan2«J>)^) 


(9a) 


y0  =  n  +  R2(Ri'1Yosin20  =F  COS0C1  -  Yp2W 
sin20  -  Yp2a2R^  2tan2$)^) 


(9b) 


z0  =  C  +  Yp*R2Ri"ltan<>  •  (9c) 

> 

For  equations  (9)  to  apply,  0  and  ^  are  subject  to  the  condition 
that: 


<  0,  < 

Inc 


sin  *(v/c) 
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(10) 


y 

where  v  is  equal  to  a  for  KP  and  £  for  XS. 

Transmitted  Waves 

The  transmitted  waves  are  due  to  the  numerous  combinations  of  waves 
refracted  inside  the  cylinder  before  they  are  transmitted(refracted) 
into  the  fluid  medium  outside  the  cylinder  in  which  the  recording  is 
being  done.  The  secondary  source  points  Q(|,n>C)  OQ  the  surface 
of  the  cylinder  are  redefined  as  follows  (see  Figure  12) : 

4  =  asin(0  +  0^)  (lla) 

q  =  YQ  -  acos(0  +  0i)  (lib) 

4  =  atan$  +  b^  (lie) 


b^  =  2a*sin(0^/2)tan$/YQsin0 


(lid) 


0.  =  7 r  -  2(sin"1(YiYosin0/(Yo*  +  a*  -  2aYocos0)*) 


(lie) 


where  i  =  1  for  P  waves  and  2  for  S  waves  and 


Yj  «  «/c;  y2  =  0/e;  6ti  *  •i“-i^Yi*i“8inc) 


(Ilf) 


For  a  given  ray  path  defined  by  an  initial  6  and  the  three 
dimensional  ray  paths  interior  to  the  elastic  cylinder  trace  out  a 
helix  on  the  surface  of  the  cylinder.  Figure  12  defines  the  parame¬ 
ters  to  be  used  to  obtain  equation  (lid).  Figures  13,  14  and  15 
show  thw  plan  views  of  the  three  triangles  used  to  obtain  the  fol¬ 
lowing  relationships  leading  to  equation  (lid). 

SQ  «  R  s  (YQ2  -  2a¥0cos9  +  a2se.c2$)*  (12) 

SC  =  Di  =  (YQ2  +  a2tan2$)*  (13) 

QC  =  a  =  cylinder  radius  (14) 


(15) 


30 


Figure  12:  Transmitted  wave  parameter  description 
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Using  equations  (13),  (14)  and  (15)  we  can  rewrite 

equation  (12)  and 

i 

solve  for  cosi|t^: 

R2^  =  D2^  +  a2  -  2aD^cosiJ(j 

(16) 

; 

t 

cos^  =  -Y0cos6(Y20  +  a2tan2$)”^ 

i 

(17) 

c 

t2  =  HC  »  cosijjj  =  a/T2 

»-» 

00 

< 

=  a(^2Q  +  a2tan24»)^/YQsin0 

f. 

D9)  [ 

< 

Tj  =  HQ  =  (T22  -  a2)** 

[ 

(20)  j: 

!’ 

c 

Tj  =  a(Y2Qsin20  +  a2tan2$)^/YgCos0 

W 

i 

A 

(21) 

t- 

T3  =  atan0 

■  " 

i- 

i 

(22)  L 

t 

tan^  =  (T2j  -  T23)%/T3 

(23) 

c 
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r 

f 

* 

-V 

(24) 


tani^  =  atan$/YgSin0  . 


bA  =  QBtaniJ»2 


(25) 


QB  =  2asin(6i/2)  . 


(26) 


Using  equations  (24),  (25)  and  (26)  we  obtain  equation  (lid)  for  the 
pitch  of  the  helix:  For  the  first  transmitted  P  wave,  KPK,  the 
expanding  secondary  source  radius  is: 


R^  =  ct  -  R^  -  2(c/cr)a*sin(01)(l+  a2tan2<J>/Yg2sin20)^  (27) 


Similarly,  the  expanding  radius  for  the  first  transmitted  S  wave,  KSK,  is 


R^  =  ct  -  R^  -  2(c/P)a2sin(02)(l+  a2tan2<f>/Yp2sin20)^  (28) 

Substituting  equations  (27)  and  (28)  into  equations  (3)  and  (4),  we 
obtain  the  parametric  wavefront  and  ray  path  equations  for  KPK  and 


KSK  as  follow: 


/ 

*0  =  i  -  ((Zj  +  Z2)/(l  +  Yi))cos(e  +  0.)  ±  sin(0  +  0.) 

<Rj2  -  ( (Zj  +  Z2)/(l  +  Yi))2  -  (2q  -  02)h  (29a) 

y0  =  H  ’  ((Zj  +  Z2)/(l  +  Y.))sia(0  +  0i)  f  cos(0  +  0.) 

(Rj2  ‘  ((Z1  +  Z2)/(1  +  Yi))2  “  (20  '  £)2)%  (29b) 

ZQ  ~  C  “  Rjatan^(Ri"1(2ac/ci)sin(0i/2)(l  +  a2taa2<ji//Y02 

sia20)  Js)/(YQsin0  +  2asin(0i/2))  (29c) 

Z1  =  Rj  C(-*c/cj)  (1  +  a2tan2<j>/Y()2)Vcos(0i/2)  +  (2a3c/c. 

J 

Yo2)sin(0i/2)taa2<Kl  +  a2taa2<j>/Yo2sia20)"^csc20cot0 
+  YoRi"lsin0)  (29d) 

z2  =  (20  -  C)(aYitaa^cos(0i/2)/Y()sia0  -  (a/YQ)taa$ 

csc0cot0sia(0^/2))  (29e) 


ladaSSSI 
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Y^  =  (YiaYQ2sin26  -  y^YgCosS^g2  +  a2  -  2aYQcos0)) 

/((Y02  +  a2  -  2aYocos0)(Yo2  +  a2  -  2aY0cos6  -  Y^V 

sindfy  .  (29f) 

For  the  KPK  wavefronts  and  ray  paths  use  equations  (28)  and(29)  with  i 
=  1  and  j  =  4.  Similarly  for  KSK  use  i  =  2  and  j  =  5  in  equations 
(28)  and  (29). 

Multiply  internally  reflected  and  converted  transmitted  waves 

To  obtain  these  wavefront  and  ray  path  parametric  equations,  we 
again  must  redefine  the  secondary  source  points  Q(£>H>C)  which 
are  on  the  surface  of  the  cylinder  as: 

£  =  asin(8  +  10^  +  mO^)  (30a) 

* 

*  / 

q  =  Yq  -  acos(0  +  10j  +  mO^)  (30b) 

£  =  atan$  +  10^  +  (30c) 
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x0  =  t  -  ((Qj  +  Q2)/(l  +  lYj  +  mY2))cos(0  +  1B1  +  m02) 

±  sin(0  +  10J  +  n02)(R1>n2  "  (z  -  O2  -  CCQj  +  Q2)/(l  +  lYj 

+  mY2)))*  (32a) 

y0  =  H  -  ((Qx  +  Q2)/0  +  lYj  +  mY2))sin(0  +  10j  +  ®02) 

+  cos(e  +  iex  +  n02)(R1)B2  -  C Zq  -  C)2  -  C(Q1  +  q2)/(i  + 

lYj  +  mY2)))* 


(32b) 


z0  =  i  +  R1>m(aRi_1Y0tan<Mne  -  2aztan0( (lc/aJsinO^ 


2))(1  +  a2tan2«)i/Y02sin2e)^)/Y0sine(Y0sine  +  23(5111(6^2)  + 
sin(82/2)))  (32c) 


Q1  =  (20  ‘  O((a2tan2<f./Yosine)(Y1cos(01/2)  +  Y2cos(02/2)) 

-  (2a2tan<|>/Yo)(sin(01/2)  +  sin(02/2)))  (32d) 

Q2  =  Ri>n.(YORi"lsin0  +  ((10/0)510(6^2)  +  (mc/p)sin(02/2)) 

(1+  a2tan2<J>/Yo2sin20)"^(2a3/Yo2)tan2(()csc2ecot0  + 

y 

(1  +  a2tan2<)i/Yo2sin2e)((lc/a)Y1cos(01/2)  +  (mc/p)Y2 

cos02/2)))  .  (32e) 

Yj  and  Y2  are  given  by  equation  (29e).  Equations  (32)  will  pro¬ 
vide  the  wavefronts  and  ray  paths  for  all  of  the  transmitted  waves. 
If  1  =  1  and  m  =  0  or  1  =  0  and  m  =  1  are  applied  to  equations  (32)  we 
obtain  equations  (17)  for  KPK  and  KSK,  respectively.  If  we  let  (j»  = 

0  in  equations  (29)  we  obtain  the  two  dimensional  equations  as  fol¬ 
lows: 
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y 

*  =  K  -  ((2a  +  Z2)/(l  +  Yi))cos(6  +  6.)  ±  sin(6  +  6.) 

(rj2  "  ((Z1  +  Z2)/(1  +  Yi))2  '  (z  "  (33a) 


y  =  n  -  ((Zj  +  z2)/(i  +  Yi))sin(e  +  ei)^r  Cos(e  +  ei) 


(*j*  -  ((Zj  +  Z2)/(l  +  Y.))2  -  (2  -  C)2)*  (33b) 


2  =  0 


(33c) 


•1  . 


Z1  =  rj(YOri  sin6  "  C*c/ci)Tico»(ei/2)) 


(33d) 


z2  =  o 


(33e) 


Y^  =  see  equation  (29f) 


(33f) 


Cj  =  ct  -  x±  -  (c/ci)2asin(6i/2) 


(33g) 


ri  *  <V  +  *2  '  2aYocos8)J|  ’ 


(33h) 
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The  following  Figure*  16  and  17  display  the  wavefronts,  ray  paths 
and  caustics  generated  by  the  various  transmitted  waves  according  to 
the  two  dimensional  equations  obtained  from  equations  (32)  by  allow¬ 
ing  $  to  become  0.  These  waves  are  identified  in  the  data  in  Fig¬ 
ures  7  through  10.  All  transmitted  converted  waves,  which  have  the 
same  number  of  P  and  S  ray  segments  interior  to  the  cylinder,  arrive 
simultaneously  at  the  receivers  ,.e.g.,  KPSK  and  KSPK  and  K2PSK, 
KPSPK  and  KS2PK  arrive  superimposed  as  single  wavelets.  Figure  18 
shows  the  ray  combination  for  the  simultaneous  arrival  of  KPSK  and 
KSPK.  The  interpreted  data  are  labeled  KIPmSK  since  these  waves 
arrive  simultaneously.  The  form  of  the  transmitted  mode  converted 
waves,  KIPmSK,  varies  greatly  between  receivers  as  shown  in  Figures 
7  through  10.  However,  the  K1PK  and  KmSK  homogeneous  transmitted 
waves  retain  their  waveform  between  receivers.  The  variance  of  the 
waveform  character  for  the  KIPmSK  events  and  the  retainment  of  the 
waveform  character  for  the  K1PK  and  KmSK  events  is  noticeable  by  an 
inspection  of  the  data  received  in  Regions  I  and  II  as  displayed  in 
Figures  7  and  8.  Transmitted  waves  are  received  in  all  four 
regions.  Their  interference  with  the  diffracted  waves  is  most  not¬ 
iceable  in  Regions  I  and  II.  It  is  within  these  two  regions  that 
the  diffractions  retain  enough  energy  so  that  they  are  not  masked  by 
the  transmitted  waves  as  they  are  in  Regions  III  and  IV. 


KPSK  ' - 

t  =  80  usees] 


d.)  KPSPK 

t  a  80  usees. 


_)  KSp<  , - - - c7TK2PSlT 

t  *  80  usees.  =  80  usees 


e.)  KSPSK 

t  *  80  usees. 


Figure  17:  Transmitted  wavefronts 
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Diffracted  Waves 

The  path  of  a  diffracted  wave  on  the  cylindrical  acoustic-ealstic 
interface  is  a  helix,  as  shown  in  Figure  19,  beginning  at  the  point 
of  excitation.  The  slope  of  the  helix  is  defined  by  the  tangent  to 
the  cylinder  at  this  excitation  point.  The  slope  is  given  parame¬ 
trically  as: 


bdk  =  a2tan$/Yflsineck  (34) 

where  ©ck  *  sin  *(c/ck),  k  =  1,2,3.  For  k  =  1  use  a,  for 

k  =  2  use  0  and  for  k  =  3  use  v  .  For  the  KK.K  wave  the  helix 
r  r  d 

slope  is : 

bdk  =  aztan$/Ric.  (35) 

The  length  of  the  incident  acoustic  ray  path  for  all  of  the  dif¬ 
fracted  waves  is: 

Ric  =  (Yq2  -  2aYosin0c  +  a2sec2<»*  (24) 

where  6c  and  satisfy  the  following  equations  for  the  dif¬ 
fracted  waves: 
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Figure  Diffracted  wave  helical  ray  path  on  the 
cylinder. 
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For  KP ,K: 
a 


7T  -  cos‘1((YQ2  -  *Y0cos6c  +  a2sec2<|>c)/aRic))  = 


sin  *(c/ci) 


For  KS .K: 
a 


7j T  -  cos*‘1((Yq2  -  aYQcos0c  +  a2sec2$c)/aRic))  = 


sin-1(c/p) 


For  KRdK: 


J(  -  cos  1((Yq2  -  aYQcos8c  +  a2sec2$c)/aRic))  = 


sin-1(c/vr) 


For  KKjK: 


cos"1^2  -  *Yocos0c  +  a2sec2*c)/aRi(;))  =  if  /2 


The  parametric  equations  describing  the  positions  of  the  secon¬ 
dary  source  points  Q(£,H»C)  on  the  cylindrical  acoustic-elas¬ 
tic  interface  are  given  by: 

£  =  asin6  (40a) 

H  =  Y0  -  acos0  (40b) 

£  =  atan$  +  bdk(0  -  0c).  (40c) 

The  expanding  radii  from  the  secondary  source  points  for  each  dif¬ 
fracted  wave-type  are  given  by: 

R(KP.K)  s  ct  -  R.  -  sc (0  -  0  )/ct  (41a) 

a  ic  c 

R(KSdK)  =  ct  -  Ric  -  sc (0  -  0c)/p  (41b) 

R(KRdK)  =  ct  -  Ric  -  sc(0  -  9C)/Vr  (41c) 

R(XKdK)  =  ct  -  Ric  -  sc (0  -  0c)  (4ld) 
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where  s  equals  the  arc  length  of  the  helix  for  1°  rotation. 


s  =  (a2  +  p2)^  (4le) 

Substituting  equations  (40)  and  (41)  into  equations  (3)  and  (4)  we 
obtain  the  Substitution  of  equation  (26)  for  p  in  equation  (4le) 
yields  the  arc  length  for  a  1°  rotation. 

s  =  a(YQ2sin2ec  +  (0. 1746)2tan2<frc)  ^/YgSine^,  (4lf) 

Substituting  equations  (40)  and  (41)  into  equations  (3)  and  (4)  we 
obtain  the  following  wavefront  and  ray  path  parametric  equations  for 
the  diffracted  waves: 


For  KP  .K  : 
a 


X0  =  i  +  AjCosO  ±  sin8(R2(KPdK)  -  (z„-  £)2  -  Aj2)*  (42a) 


Y0  ■  H  +  AjSinO  +  cos0(R2(XPdK)  -  (z0-  C)2  -  A^)5* 


(42b) 


(42c) 


zQ  =  C  +  RCKP^JCaaR.^^Qtancfisine  +  0  -  0c)/a(YQ 

\ 

sinG  +  0  -  0  ) 
c  c 


Aj  =  (c/a)R(KPdK)  -  (z,-  C)tan$/Y0sinec 


(42d) 


Similarly  to  obtain  the  equations  necessary  to  calculate  the 
wavefronts  and  ray  paths  for  XS^X,  KR^K  and  XX^X  we  need  only 
to  substitute  equations  (41),  along  with  the  correct  velocity  in 
place  of  a,  into  equations  (42). 

To  obtain  the  two  dimensional  equations  for  the  diffracted  waves 
we  set  $  ~  0  in  equations  (33-42).  Doing  so  we  obtain: 


4  =  asinfl  (43a) 

H  =  YQ  -  acos0  (43b) 

C  =  0  (43c) 

ric  =  <V  +  **  •2aYoCO80c)^  (43d) 
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(44a) 


y 

Xq  =  5  +  (rc/a)cos8  ±  rsind(l  -  c2/ct2)^ 

7q  =  H  +  (rc/a)sin8  ^  rcos8(l  -  c2/a2)3j  (44b) 

z0  =  0  (44c) 

r  =  r(KP<jK)  =  ct  -  ric  -  (ac/a)(8  -  8c)  (44d) 

By  substituting  the  appropriate  velocity  in  place  of  a  in  equa¬ 
tions  (33)  they  can  be  transformed  into  the  equations  describing  the 
wavefronts  and  ray  paths  for  the  remaining  diffracted  waves.  The 
diffracted  wavefronts  for  all  four  diffracted  waves  are  shown  in 
Figure  20  as  computed  from  the  above  equations  at  t  =  70  psec.. 

Analysis 

The  events  labeled  in  the  data,  Figures  7,  8,  9  and  10,  were  picked 
by  calculating  the  arrival  times  of  the  different  wavefronts  at  the 
receivers  using  the  simple  ray  theory  discussed  in  equations  1  through 
44  of  the  previous  section.  Since  the  velocities  of  the  acoustic  and 
elastic  media  were  known,  this  information  was  also  used  in  the 
interpretation  of  the  arrivals  of  the  diffracted  waves.  By  placing 
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the  receivers  on  concentric  arcs  about  the  cylinder ,  the  identifica¬ 
tion  of  the  diffracted  waves  is  simplified  because,  for  a  given  type 
of  diffracted  wave  received  at  two  receivers  on  the  same  arc,  the  only 
part  of  the  ray  path  that  changes  is  the  path  on  the  cylinder.  Only 
this  segment  of  the  entire  ray  path  changes  because  the  diffracted 
wave  is  always  radiated  from  the  cylinder  at  the  same  angle.  There¬ 
fore,  the  time  interval  between  the  arrival  of  the  diffracted  wave  at 
two  receivers,  0^  and  0^,  on  the  same  arc  can  be  calculated  in 
the  following  manner: 

AT  =  sCOj  -  62)/V  (45) 

where  s  is  the  helical  arc  on  the  cylinder  and  V  is  the  velocity  of 
the  diffracted  wave  on  the  cylinder.  If  the  arrival  times  of  the  dif¬ 
fracted  events  are  plotted  as  a  function  of  0,  then  they  will  lie  on 
straight  lines  whose  slopes  will  be  proportional  to  the  velocity  of 
the  diffracted  wave.  This  phenomenon  is  most  easily  seen  in  Regions  I 
and  II  as  shown  in'  Figures  7  and  8  for  both  2D(coplanar)  and 
3D(non-coplanar)  data.  For  non-symraetric  receiver  positions  these 
diffracted  events  split  into  front  and  back  diffractions  propagating 
in  both  directions  around  the  cylinder.  The  time  interval  between 
arrival  at  the  same  receiver  is  given  by 

^ST  =  2s/v,  (46) 
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where  s  is  the  arc  length  of  the  helical  path  on  the  cylinder  and 
v  is  the  velocity  of  the  diffraction.  Again  this  phenomenon  is  best 
seen  in  the  dats  collected  in  Regions  I  and  II  (see  Figures  7 
and  8)  . 

Interpretation  and  Discussion 

The  following  events  are  identified  in  this  Region  I: 

1. )  KPK,  the  transmitted  P  wave 

2. )  KP-K*3,  the  back  diffracted  P  wave 

d 

3. )  KSK,  the  transmitted  S  wave 

f  b 

4. )  KS.K  '  ,  the  front  and  back  diffracted  S  waves 

d 

5. )  KRdKf'b,  the  front  and  back  diffracted  Rayleigh  waves 

6. )  K3PK,  the  twice  internally  reflected  P  wave 

f ,  b 

7. )  KK.K  ,  the  front  and  back  diffracted  acoustic  waves 

a 


8. )  the  first  circumnavigated  front  and  back 

diffracted  P  waves 

9. )  K5PK,  the  four  time  internally  reflected  P  wave 

10. )  KPPSK,  KPSPK,  KSPPK,  labeled  K2PSK 

11. )  KPSSK,  KSPSK,  KSSPK,  labeled  KPSSK 

The  data  for  Region  I,  the  shadow  zone,  are  displayed  in  Figure  7 
For  this  discussion  the  coplanar  case  will  referred  to  as  2D  and  the 
non-coplanar  case  as  3D.  Both  the  2D  and  the  3D  data  exhibit  the  same 
phenomena.  The  KPK  event  is  the  first  arrival  followed  closely  by  the 
KP^K  event.  At  the  symmetrical  receiver  position  0°  azimuth,  the 
KPK  event  is  clearly  identifiable.  Since  the  receivers  are  placed  on 
a  concentric  arc  about  the  cylinder  the  moveout  of  the  diffracted  wave 
is  constant(see  equation  (45))  and  the  arrivals  should  appear  to  lie 
on  a  straight  line  when  the  arrival  time  is  plotted  as  a  function  of 
azimuth.  This  is  easily  observed  in  both  the  2D  and  the  3D  data  (see 
Figure  7).  The  KP^K^back  diffracted  P  wave)ev«nt  travels  a  lon¬ 
ger  path  to  arrive  at  receivers  5°  to  30°  and  ,  therefore,  arrives 

later  as  the  receiver  azimuth  increases.  The  KP,K*(front  dif- 

o 

fracted  P  wave)  event  travels  an  increasingly  shorter  path  to  these 

same  receivers.  However,  due  to  the  magnitude  of  the  KPK  event  the 

KP^K*  event  can  not  be  detected.  Due  to  this  splitting,  the 

KP.K^’k  events  destructively  interfere  at  azimuth  5°  for  both 
0 

cases  and  at  10°  for  the  2D  case.  The  next  strong  arrival  is  the  KSK 
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event.  This  event  is  very  prominent  but  as  the  receiver  nears  the 
shadow  zone  boundary,  the  KSK  and  the  KP^K^  events  begin  to  inter- 


fere.  At 

25° 

the  KSK  event  is 

lost 

but  the  diffracted  event 

is  iden- 

tifiable. 

At 

30°  both  events 

are 

identifiable.  The  front 

and  back 

traveling  diffractions  are  most  observable  by  tracking  KS^K^’^  and 

KR.K*’k.  These  events  form  a  V  as  the  receiver  azimuth  increases, 
d 

Interference  between  the  KSK  and  the  KS^K^  events  begins  to  occur 
when  the  source  configuration  is  no  longer  symmetric.  It  becomes 
increasingly  difficult  to  detect  the  correct  arrival  time  of  KS^K^ 
and  at  25°  this  interference  destroys  both  events  creating  a  data 
void.  K3PK  is  easily  detectable,  but  as  the  shadow  zone  boundary  is 
approached,  interference  with  KR^K(the  back  diffracted  Rayleigh 
wave)  increases.  This  interference  is  so  great  that  it  appears  that 
only  one  event  arrives.  As  the  shadow- boundary (30°  azimuth)  is 
crossed  ,  the  KR^K^  and  K3PK,  as  well  as,  the  KP^jK^  and  KSK 
events  cross  and  become  separately  identifiable  again.  The  multiply 
internally  reflected  mode  converted  transmitted  waves  K2PSK  are  evi¬ 
dent  in  the  time  interval  between  K3PK  and  KK^K.  Again,  the  dif¬ 
fracted  acoustic  wave,  KK^K^,  crosses  K2PSK  at  azimuth  20°  and  it 
also  crosses  K3PK  and  KR.K^5  at  30°  azimuth(  the  shadow  zone  boun- 

Q 

dary).  Beyond  30°  azimuth,  KK^K^  is  not  detectable  until  it  cir- 
cum-navigates  the  cylinder.  Diffractions  which  circum-navigate  the 
cylinder  are  detectable  as  evidenced  by  the  arrival  of  KP^K1  fol¬ 
lowing  the  KK^K  event  at  0°  azimuth.  The  later  arrivals  consist  of 
KP2SK  and  K5PK.  Other  KIPmSK,  K1PK  and  KmSK,  as  well  as  other  circum- 
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navigating  diffractions,  are  extremely  difficult  to  interpret  with  any 
confidence  so  they  are  not. 

The  following  events  are  identified  in  this  Region  II: 

1. )  KPK,  the  transmitted  P  wave 

2. )  KP^K^,  the  back  diffracted  P  wave 

3. )  KSK,  the  transmitted  S  wave 

4. )  KS^K^’k,  the  front  and  back  diffracted  S  waves 

5. )  KR^K^’k,  the  front  and  back  diffracted  Rayleigh  waves 

6. )  K3PK,  the  twice  internally  reflected  P  wave 

7. )  KK^K^’k,  the  front  and  back  diffracted  acoustic  waves 

8. )  KP^K1^’^,  the  first  circumnavigated  front  and  back 

diffracted  P  waves 

9. )  K5PK,  the  four  time  internally  reflected  P  wave 

10. )  KPPSK,  KPSPK,  KSPPK,  labeled  K2PSK 

11. )  KPSSK,  KSPSK,  KSSPK,  labeled  KPSSK 

In  Region  II  the  first  phenomenon  noticed  is  the  amplitude  decay  of 
KPK  and  KSK.  This  is  not  as  noticeable  in  Region  I.  Each  of  the 
other  events  which  are  identified  for  the  shadow  zone  are  also 
observed  and  identified  in  this  region.  Besides  these  events,  Region 
II  adds  KFSK,  KSPK  and  K  to  this  list.  The  first  arrival  in  Region  II 
is  the  KPK  event.  Although  these  reveivers  are  in  the  illuminated 
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zone,  diffracted  P,  S  and  Rayleigh  waves  are  still  received.  The 
KP^K^  event  is  not  evident  since  it  arrives  in  the  time  window  of 
KPK.  The  KP^K^  arrival  at  35°  is  superimposed  on  KSK  and  remains 
unseparated  from  it  until  azimuth  50°.  For  receivers  between  at  azi¬ 
muths  35°  and  50° ,  four  different  events  are  superimposed  within  the 
time  window  98.25  to  103.25  psecs.,  making  identification  of  the 
separate  events  difficult.  These  events  are  KPdKb,  KSK,  KSdKf 
and  KRdK^.  The  KSdK^  event  is  not  visible  as  a  separate  event 

at  35°  azimuth  but  KR.K^  has  an  identifiable  arrival  time  but  no 

a 

termination  point.  Arriving  on  the  tail  of  KRdK*  is  the  mode  con¬ 
verted  events  KPSK.  These  mode  converted  arrivals  at.*  unidentified  in 
the  shadow  zone.  Another  event  which  is  non-existent  in  the  shadow 
zone  ,  but  does  exist  in  the  near  shadow  zone  and  in  the  remainder  of 
the  illuminated  zone,  is  the  direct  acoustic  wave,  K.  K  arrives  at 
receiver  35°  at  112  psecs.  coupled  with  KRdKb.  At  40°,  the  K 
event  arrives  coupled  with  KSdKb  and  K3PK.  Since  the  radiation 
pattern  of  the  source  transducer  is  narrow(Figure  3),  the  K  event  is 
low  in  amplitude  and  shows  no  evidence  of  its  existence  at  45°  ?od 
50°.  At  55°,  the  K  event  arrives  coupled  with  the  mode  converted 
events  KPSK.  Throughout  the  near  shadow  region  the  diffracted  events 

KS.Kb  and  KR.Kb  are  identified.  Mode  converted  events  K2PSK 
a  fl 

and  KP2SK  are  detected  and  identified  in  this  region.  The  KP^jK1^ 
arrival  in  this  region  displays  a  moveout  to  earlier  times  as  the 
receiver  azimuth  increases.  It  should  be  noted  that  back  diffracted 
acoustic  waves  are  present  in  the  illuminated  zone,  as  well  as,  the 
diffracted  waves  KPdKf,b,  KSdKf,bK  and  KRdKf,b. 
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The  following  events  are  identified  in  Region  III: 


1. )  K,  the  direct  acoustic  wave 

2. )  KPK,  the  transmitted  P  wave 

3. )  KP^K^’k,  the  front  and  back  diffracted  P  waves 

4. )  K2PK,  the  first  internally  reflected  P  wave 

5. )  KSK,  the  transmitted  S  wave 

6. )  KS^K^'k,  the  front  and  back  diffracted  S  waves 

7. )  KR^K^’k,  the  front  and  back  diffracted  Rayleigh  waves 

8. )  KK  ,  the  reflected  acoustic  wave 

9. )  K3PK,  the  second  internally  reflected  P  wave 

10.)  KPSK,  KSPK,  labeled  KPSK 

12. )  K2SK,  the  first  internally  reflected  S  wave 

13. )  KPjK1^,  the  first  circumnavigated  diffracted  P  wave 

14. )  KPSSK,  KSPSK,  KSSPK,  labeled  KPSSK 

16. )  K3SK,  the  second  internally  reflected  S  wave 

17. )  KPSSK,  KSPSK,  KSSPK,  labeled  KPSSK 

18. )  KPSSK,  the  mode  converted  P  to  SS  wave 

19. )  KPSSSK ,  KSPSSK,  KSSPSK,  KSSSP,  labeled  KP3SK 

In  this  region  certain  events  are  received  and  identified  which  are 
not  present  in  the  shadow  or  the  near  shadow  regions.  Since  receiver 
90°  and  105°  are  within  the  reflection  region,  reflected  acoustic 
waves,  KKr,  are  received  and  since  receiver  75°  is  outside  of  this 
region  no  KKf  event  exists(see  Figure  3).  Transmitted  waves  with  an 
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even  number  of  P  ray  paths  and  S  ray  paths  without  any  mode 
conversions  are  observed  at  these  receivers.  K2PK  and  K2SK  are  uni¬ 
dentified  at  receivers  in  the  shadow  and  near  shadow  zones  but  they 
are  observed  for  receivers  at  75°,  90°  and  105°  azimuth.  A  splitting 

effect  for  K2PK  and  K2SK  is  evident  in  this  region  caused  by  internal 
reflections  traveling  in  both  directions  around  the  cylinder.  On  exa¬ 
mination  of  the  transmitted  ray  diagrams  in  Figures  16  and  17,  this 
effect  is  seen  if  it  is  remembered  that  these  diagrams  plot  only  one 
half  of  the  ray  path  and  wavefront.  Dual  arrivals  of  K2PK,  K2SK  and 
K3PK  are  identified  in  Figure  9.  Host  of  the  observed  events  in  this 
region  are  combinations  of  events.  The  first  arrival  is  a  combination 

of  K,  KPK,  KP^K^  and  K2PK  at  receiver  75°.  At  90°,  the  K  event 

arrives  much  earlier  than  the  combination  KPK,  KP,Kf  and  K2PK  with 

a 

such  a  low  amplitude  that  it  is  unrecorded.  At  all  three  receivers, 

KSK,  KS^K*  and  KR^K^  arrive  as  overlapping  wavelets.  From 

trace  to  trace  the  form  of  the  signal  changes  because  each  of  these 
events  has  a  different  moveout  such  that  they  cross  each  other.  At 
105°,  both  diffracted  events,  KS^K^  and  KR^K*,  arrive  before 
KSK  due  to  their  larger  moveouts  as  the  receiver  azimuth  increases. 
Between  116  and  120.5  psecs.,  the  K2PK,  the  K3PK  and  the  KPdKb 
events  arrive  very  close  together.  At  receiver  azimuth  90°,  the 
KPdKb  event  is  undetected  but  is  reappears  at  105°.  K3PK  arrives 
at  receiver  90°  on  the  tail  of  the  reflected  acoustic,  KK  .  K2PK  is 
present  on  all  three  traces.  Within  the  time  window  123  to  128.5 
psecs.  at  receiver  75°,  the  KPSK  events  arrive  simultaneously  fol- 
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lowed  by  a  second  K3PK  event  due  to  the  splitting  of  transmitted 
waves,  previously  mentioned.  At  90°,  KPSK  and  K3PK  arrive  closer 
together  and  at  105°  they  combine  to  form  a  single  low  amplitude 
event.  K2SK  arrives  following  this  group  at  all  three  receivers  fol¬ 
lowed  by  the  mode  converted  events  K2PSK.  A  second  K2SK  event  is 
identified  at  receiver  90°  at  134.5  psecs.  coupled  with  the  first 
circumnavigated  front  diffracted  P  wave,  KP^K**.  The  diffracted  S 
wave,  KS^K**,  appears  at  all  three  receivers  and  it  displays  its 
characteristic  moveout  arriving  later  as  the  receiver  azimuth 
increases.  The  back  diffracted  Rayleigh  wave,  KR^K^,  is  unidenti¬ 
fiable  at  75°,  but  it  is  identified  at  90°  and  105°.  It  also  displays 
moveout.  Finally,  KP2r>"  events  are  observed  arriving  at  receiver  75° 
and  KP3SK  are  observed  at  90°  arriving  simultaneously.  In  this 
region,  the  central  part  of  the  illumina'ted  zone,  more  transmitted 
waves  are  observed  than  in  the  shadow  zone  creating  a  more  complicated 
signal.  The  diffracted  waves  are  still  observed  in  this  region. 

Figures  10a  and  10b  displays  the  data  for  Region  IV  for  which  the 
source  has  been  moved  to  a  distance  of  6Rq  from  the  cylinder  and  the 
receivers  are  placed  on  a  circular  arc  of  the  same  radius  about  the 
cylinder  at  azimuths  120°,  135°,  150°  and  160°  as  shown  in  Figure  6. 

This  configuration  covers  the  deep  illuminated  zone  where  backscatter- 
ing  predominates.  Reflected  acoustic  waves,  KK^K  are  the  first  and 
most  predominant  of  the  received  events  in  this  region.  The  waveform 
is  clipped  because  its  amplitude  is  much  greater  than  the  amplitude  of 
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the  remainder  of  the  scatterred  and  diffracted  signal.  Very  little 
diffracted  energy  is  observed  in  this  region.  The  front  diffracted 
waves,  KP^K^,  KS^K^  and  KR^K^,  arrive  within  the  time  win¬ 
dow  of  the  KK^  event  and  are  therefore  undetectable  due  to  the  large 
amplitude  of  the  reflected  acoustic  wave,  KK^.  The  back  diffracted 
waves  are  undetected  because  they  are  attenuated  rapidly  and  they 
arrive  in  the  same  time  windows  as  the  larger  amplitude  transmitted 
waves . 


In  Figures  10,  the  following  transmitted  waves  are  identified  in 
Region  IV: 

1. )  K2PK 

2. )  KPSK  and  KSPK,  labeled  KPSK  - 

3. )  K3PK 

4. )  KPSPK,  KPPSK  and  KSPPK,  labeled  K2PSK 

5. )  KSSK 

6. )  KPPPSK,  KPPSPK ,  KPSPPK  and  KSPPPK,  labeled  K3PSK 

7. )  K4PK 

8. )  K3SK 

9. )  KSPSPK,  KPSSPK , KSSPPK , KPPSSK  and  KSPPSK,  labeled  K2P2SK 

10.)  KSSPK,  KSPSK  and  KPSSK,  labeled  KP2SK 


Just  as  the  diffracted  waves  split  into  two  arrivals  propagating  in 
both  directions  around  the  cylinder,  a  similar  splitting  occurs  for 
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the  transmitted  waves.  The  K3PK  event  shows  a  dual  arrival  at  160° 


overlapping  each  other  at  190.5  psecs..  The  mode  converted  waves 
K3PK  display  this  phenomenon  on  all  traces  with  arrivals  at  210.8  and 
217.6  psecs.  at  receiver  120°.  These  dual  arrivals  become  closer 
together  as  the  receiver  azimuth  increases  as  seen  in  Figures  10. 
This  backscattered  region  contains  mostly  transmitted  waves  . 


I 


% 


l 


The  amplitude  spectra  of  the  data  have  many  peaks  and  troughs 
unlike  the  source  wavelet(see  Figure  2b).  Modulation  of  the  amplitude 
spectrum  is  caused  by  the  multitude  of  arriving  events  within  the  time 
window  analyzed.  The  peaks  are  related  to  the  relative  arrival  times 
between  two  events  The  reciprocals  of  the  peak  frequencies  are  the 
relative  arrival  time  difference.  Having  already  interpreted  the 
data,  comparison  of  the  spectral  peaks  to  the  relative  arrivals  of  the 
data  is  done.  Tables  I  and  II  show  the  correspondence  between  the 
spectral  peaks  and  the  relative  relative  arrival  times  between  events 
for  the  2d  coplanar  case  at  receivers  0°  and  5°  in  the  shadow  zone. 
Figure  29  displays  the  spectra  for  receivers  0°  and  5°  for  the  data 
presented  in  Tables  I  and  II. 

CONCLUSIONS 


A  three  dimensional  model  experiment  has  been  performed  to  study 
the  acoustic-elastic  wave  scattering  and  diffraction  of  a  transient 
acoustic  wave  emitted  from  a  point  source  by  an  infinite  elastic  cyl¬ 
inder  embedded  in  a  fluid  medium.  This  is  the  first  three  dimensional 
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model  experiment  using  multiple  receiver  positions  throughout  all  of 
the  scattering  regions.  The  data  verify  the  existence  of  KPK,  KSK, 
K2PK,  K2SK,  KPSK ,  KP2SK,  K2PSK,  K3PK,  K3SK,  K3PSK,  KP3SK,  K5PK,  K, 
KKr,  KPdf,bK,  KPdlfK,  KSdf,bK,  KRdf>bK,  KKdf,bK 

and  those  transmitted  converted  waves  which  have  the  same  number  of  P 

and  S  wave  ray  segments  interior  to  the  cylinder  as  those  stated  here. 

The  verification  of  the  existence  of  KP.^,bK,  KS,^'bK, 

a  d 

KRd^’bK,  KKd^’bK,  KPd^K  in  three  dimensions  and  of 
KPd^,bK  and  KSd^,bK  in  two  dimensions  is  accomplished  for  the 
first  time.  It  is  also  for  the  first  time  that  all  of  the  diffrac¬ 
tions  and  all  of  the  transmitted  waves  have  been  shown  together  on  the 
same  data  trace  so  that  their  interaction  could  be  discerned.  By 
establishing  the  existence  of  the  various  transmitted  and  diffracted 
waves  for  a  non-coplanar(three  dimensional)  source-receiver  configura¬ 
tion,  the  verification  of  the  helical  path  traced  on  the  cylindrical 
surface  by  the  various  as  predicted  by  theory  is  accomplished.  It  is 
also  shown  that  the  diffractions  due  propagate  into  the  shadow  zones 
due  to  the  fact  that  they  propagate  in  all  directions  around  the  the 
cylinder.  Through  spectral  analysis  the  spectral  maxima  are  shown  to 
correspond  to  the  relative  arrival  time  separartion  between  various 
events.  The  non-abrupt  nature  of  the  acoustical  shadow  boundary  is 
shown  by  use  of  the  complex  cepstrum.  The  maximum  phase  component  of 
the  cepstrum  becomes  greatly  changes  in  the  vicinity  of  the  geometri¬ 
cal  shadow  boundary.  If  the  geometrical  shadow  zone  is  being  entered 
then  the  negative  cepstrum  becomes  more  complicated  and  the  terms 
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increase  in  amplitude.  However,  upon  leaving  the  geometrical  shadow 
the  negative  cepstrum  decreases  in  complexity  and  the  amplitude  of  the 
terms  decrease  also.  Envelope  detection  is  employed  to  determine  the 
arrival  times  and  the  duration  and  the  duration  of  the  arriving 
events.  From  the  envelope  detection  method  the  contamination  of  arri¬ 
vals  due  to  overlapping  is  readily  seen.  Homomorphic  filtering  is 
used  to  attempt  to  remove  the  wavelet  effects  and  decouple  the  over¬ 
lapping  events,  but  it  does  not  perform  satisfactorily  since  the  con¬ 
tamination  is  by  more  than  one  or  even  two  events  in  most  instances. 


y 
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